利用 Sympy 符号积分计算向量函数的 \(L^2\) 投影

问题表述

令四面体单元 \(K\) 的四个顶点坐标为\[ (0,0,0),\ (1,0,0),\ (0,1,0),\ (0,0,1),\]最低阶的 Nedelec 函数空间定义如下\[ \mathcal{R}_1(K)=\left\{ \begin{pmatrix} 1-x_3-x_2 \\ x_1 \\ x_1 \end{pmatrix}, \begin{pmatrix} x_2\\ 1-x_3-x_1 \\ x_2 \end{pmatrix}, \begin{pmatrix} x_3\\ x_3\\ 1-x_2-x_1 \end{pmatrix}, \begin{pmatrix} -x_2\\ x_1\\ 0 \end{pmatrix}, \begin{pmatrix} 0\\ -x_3\\ x_2 \end{pmatrix}, \begin{pmatrix} x_3\\ 0\\ -x_1 \end{pmatrix} \right\}\]计算函数 \(f(x):\mathbb{R}^3\to \mathbb{R}\) 在 \(K\) 上的 \(L^2\) 投影, 其中\[ f(x):=\begin{bmatrix} x_1^2x_2^2x_3^2(1-x_1-x_2-x_3)^2 \\ 0 \\ 0 \end{bmatrix}. \tag{1}\]

分析

设 \(\phi_i\) 为函数空间 \(\mathcal{R}_1(K)\) 的基函数. 那么函数 \(f(x)\) 的投影的定义可以表述如下\[ \int_{K} (Q_hf) \phi_i dx = \int_{K}f \phi_i dx, \quad \forall i=1,2,\ldots ,6. \tag{2}\]如果假定 \(Q_hf = \sum_{j=1}^6 c_j \phi_j\), 我们可以由上面的定义得到下面的线性方程组\[ A C = F, \tag{3}\]其中\[ A_{ij} = \int_{K}\phi_j \phi_i dx,\quad F_{i} = \int_{K}f \phi_i dx,\quad C_i = c_i.\]

因此, 我们需要计算系数矩阵 \(A\) 和右端项 \(F\). 分步说明如下.

(1) 导入必要的 sympy 对象.

from sympy import symbols, integrate, Matrix, zeros

(2) 利用矩阵对象定义函数空间和函数 \(f(x)\).

x1, x2, x3 = symbols('x1 x2 x3')
Base = Matrix([[1-x3-x2, x1, x1],[x2, 1-x3-x1, x2], \
    [x3,x3,1-x2-x1],[-x2,x1,0],[0,-x3,x2],[x3,0,-x1] ]).transpose()
f = Matrix([x1**2*x2**2*x3**2*(1-x1-x2-x3)**2,0,0])

(3) 定义三维积分函数.

def intK(f):
    return integrate(f, (x3,1-x1-x2), (x2, 0,1-x1), (x1,0,1))

(4) 积分得到 \(A, F\) 并求解. 注意, 其中涉及的矩阵乘法是标准的矩阵乘法, 并非 numpy 中的逐元素相乘.

M, RHS = zeros(6,6), zeros(6,1)
for i in range(6):
    for j in range(6):
        M[i,j] = Base[:,i].transpose() * Base[:,j]
    RHS[i] = Base[:,i].transpose() * f

A = M.applyfunc(intK)
F = RHS.applyfunc(intK)
C = A.solve(F)

通过计算可得 \[ C = \left[\begin{matrix}\frac{1}{415800} & 0 & 0 & – \frac{1}{415800} & 0 & \frac{1}{415800}\end{matrix}\right]^T\]从而, 利用投影的表达式可得\[ (Q_hf)(x) = \sum_{j=1}^6 c_j \phi_j = \left[\begin{matrix}\frac{1}{415800}\\ 0\\ 0\end{matrix}\right].\]

参考

[1] SymPy Tutorial.
[2] Python 科学计算: Sympy.

附录: 代码实现

# compute L2 projection of f(x) on the reference element
from sympy import symbols, integrate, Matrix, zeros

x1, x2, x3 = symbols('x1 x2 x3')
Base = Matrix([[1-x3-x2, x1, x1],[x2, 1-x3-x1, x2], \
    [x3,x3,1-x2-x1],[-x2,x1,0],[0,-x3,x2],[x3,0,-x1] ]).transpose()
f = Matrix([x1**2*x2**2*x3**2*(1-x1-x2-x3)**2,0,0])

M, RHS = zeros(6,6), zeros(6,1)
for i in range(6):
    for j in range(6):
        M[i,j] = Base[:,i].transpose() * Base[:,j]
    RHS[i] = Base[:,i].transpose() * f

A = M.applyfunc(intK)
F = RHS.applyfunc(intK)
C = A.solve(F)