椭圆方程的最小二乘三角单元Legendre Galerkin数值积分法

Least Squares Triangular Legendre Galerkin with Numerical Integral Method for Elliptic Equations

摘要: 研究三角形域上椭圆型方程的最小二乘三角单元Legendre Galerkin数值积分法。该方法基于最小二乘原理，在离散时采用Legendre Galerkin数值积分处理，使得格式既有Legendre Galerkin数值积分实施方便的优势，同时又有对应的代数方程的系数矩阵具有的对称正定的特点。对变系数部分的计算，则用Legendre-Gauss-Lobatto (或Chebyshev-Gauss-Lobatto)配置点插值处理。给出数值算例验证格式的谱精度。

Abstract: The least squares triangular Legendre Galerkin with numerical integration method for elliptic equations is investigated. The method is based on least squares formulation, but Legendre Galerkin with numerical integration is applied to the discretization in time, which combines the advantages of the Legendre Galerkin with numerical integration with the symmetric and positive definite systems being obtained by our scheme. For the calculation of variable coefficients, the Legendre-Gauss-Lobatto (or Chebyshev-Gauss-Lobatto (CGL)) collocation points can be used. Some numerical examples are given to verify the spectral accuracy and effectiveness of the proposed method.

1. 引言

本文中考虑三角形区域上如下椭圆型方程

$\{\begin{array}{l}-\Delta p+b\cdot \nabla p+cp=f\left(x,y\right),\text{}\left(x,y\right)\in {\Omega }_{\Delta },\\ p=0,\left(x,y\right)\in {\Gamma }_D,\\ n\cdot \nabla p=0,\left(x,y\right)\in {\Gamma }_N,\end{array}$ (1.1)

其中 ${\Omega }_{\Delta }$ 为三角形区域， $b={\left(b_1{b}_2\right)}^{\text{T}}$ ，以及 $b_i,c\left(i=1,2\right)$ ， $f\left(x,y\right)$ 均为给定的光滑函数，n表示沿着边界的单位法向量。

近年来，最小二乘有限元由于其结合了混合有限元法的优点已经被应用于流体力学领域中，如Stokes方程和Navier-Stokes方程等的数值求解 [1] [2] 。该方法的优点是可以避免Ladyzhenskaya-Babuška Brezzi (LBB)稳定性条件，所以等阶的插值多项式可用于所有变量，且所离散得到代数方程的系数矩阵具有对称正定的特点 [3] 。在求解时便于利用合适的迭代方法，如共轭梯度法等。为了得到高阶精度方法，有部分学者发展偏微分方程的最小二乘谱方法。文 [4] 中提出了椭圆型边值问题的拟谱最小二乘法，但是对于Chebyshev情况误差估计受其权函数影响而得不到最优。文 [5] 则研究耦合Legendre与Chebyshev方法且数值分析中避免了Chebyshev权函数的影响从而得到最优的误差估计。最小二乘拟谱方法还被应用于流体力学领域中Stokes方程的数值求解 [6] 。

谱方法由于其对于光滑解的偏微分方程有高阶谱精度，因此它与有限元法和有限差分法等一起成为偏微分方程重要的数值解法 [2] 。近年来，谱方法和拟谱方法被广泛应用于各领域数理模型的计算，如Navier-Stokes方程 [7] 和Fokker-Planck方程 [8] 等等。为了克服经典谱方法对矩形区域的依赖，基于Duffy’s三角形映射，文 [9] [10] 中发展了三角单元谱方法。文 [11] 中通过引进新的三角形映射研究了问题(4.4)的三角谱方法。该映射提供均匀的网格点映照 [12] [13] ，基于该映射，在 [13] 中考虑了三角单元谱元法，并在 [12] 中研究了相应的逼近理论与构造椭圆方法的三角谱元法。近年来，文 [14] 研究了耦合三角单元和矩形元谱元法。文 [15] 发展了椭圆方程和Stokes方程的三角谱元法且推广到三棱柱谱方法 [16] 。本文中利用文 [11] 的三角形映射，我们研究椭圆方程(4.4)的最小二乘三角单元Legendre Galerkin数值积分(LST-LGNI)解法。

2. 记号与三角形映射

对于任意非负整数s，设 $H^s\left(\Omega \right)$ 为经典的Sobolev空间，其中 ${\Vert \cdot \Vert }_s$ 和 $\left|\cdot \right|$ 分别表示范数和半范数。当 $s=0$ 时，则 $L^2\left(\Omega \right)=H^s\left(\Omega \right)$ ，且 ${\left(\cdot ,\cdot \right)}_{\Omega }$ 和 ${\Vert \cdot \Vert }_{\Omega }$ 分别为它的内积与范数 $L^2\left(\Omega \right)$ 。为了简单起见，记

${\Omega }_{\Delta }:=\left\{\left(x,y\right):0<x,y,x+y<1\right\},{\Omega }_{\square }:=\left\{\left(\hat{x},\hat{y}\right):-1<\hat{x},\hat{y}<1\right\}.$

对于任意的 $\left(x,y\right)\in {\Omega }_{\Delta }$ ，文 [11] 中给出了 ${\Omega }_{\Delta }$ 到 ${\Omega }_{\square }$ 的一个新型映射T：

$x=\frac{1}{8}\left(1+\hat{x}\right)\left(3-\hat{y}\right),y=\frac{1}{8}\left(3-\hat{x}\right)\left(1+\hat{y}\right),\forall \left(\hat{x},\hat{y}\right)\in \Omega .$ (2.1)

令 ${\omega }_0=\sqrt{{\left(x-y\right)}^2+4\left(1-x-y\right)}$ ，则T的逆映射 $T^{-1}$ ：

$\hat{x}=1+x-y-{\omega }_0,\hat{y}=1+x-y-{\omega }_0.$

如文 [11] 中得到Jcobian为： $J:=\det \left(\partial \left(x,y\right)/\partial \left(\hat{x},\hat{y}\right)\right)=\left(2-\hat{x}-\hat{y}\right)/16={\omega }_0/8$ 。因此，可定义内积为：

${\left(u,v\right)}_{{\Omega }_{\Delta }}={\left(J\hat{u},\hat{v}\right)}_{{\Omega }_{\square }}$ 或者 ${\left(\hat{u},\hat{v}\right)}_{\Omega }={\left({\omega }_0^{-1}u,v\right)}_{{\Omega }_{\Delta }}/8.$ (2.2)

文 [11] [13] 指出：由(2.2)得到若 $\frac{\partial u}{\partial n}|{\Gamma }_N$ 在 ${\Omega }_{\Delta }$ 的斜边中点 $\left(\frac{1}{2},\frac{1}{2}\right)$ 存在，则经过变换后有

${\left({\partial }_{\hat{x}}\hat{v}+{\partial }_{\hat{y}}\hat{v}\right)|}_{\left(1,1\right)}=0,$ (2.3)

设 ${\hat{\mathbb{Q}}}_N={\left[{\mathbb{P}}_N\right]}^2$ 其中 ${\mathbb{P}}_N$ 表示 $\left[-1,1\right]$ 次数不超过N多项式空间。因此，有限维空间可定义为：

${\hat{V}}_N:=\left\{\hat{v}\in {\mathbb{Q}}_N:{\left({\partial }_{\hat{x}}\hat{v}+{\partial }_{\hat{y}}\hat{v}\right)|}_{\left(1,1\right)}=0\right\}.$ (2.4)

设 ${\gamma }_N:=\left\{\left(i,j\right):0\le i,j\le N,\left(i,j\right)\ne \left(N,N\right)\right\}$ 和记 ${\Omega }_{\Delta }$ 上与 ${\hat{V}}_N$ 对应的有限维空间为 $V_N$ 。文 [11] [13] 定义拟插值 $I_N^L:C\left({\bar{\Omega }}_{\Delta }\right)\to V_N$ 为

$u_N:={\displaystyle {\sum }_{i,j=0}^N{u}_{ij}{\varphi }_{ij}\left(x_i,y_j\right)},\left(i,j\right)\in {\gamma }_N,$ (2.5)

其中 $u_{ij}=u\left(x_i,y_j\right)=\hat{u}\left({\hat{x}}_i,{\hat{y}}_j\right)$ ，但 ${\hat{u}}_{NN}$ 则由 ${\left({\partial }_{\hat{x}}{\hat{u}}_N+{\partial }_{\hat{y}}{\hat{u}}_N\right)|}_{\left(1,1\right)}=0$ 所确定和 ${\left\{\left(x_i,y_j\right)\right\}}_{i,j=0}^N$ 是 ${\Omega }_{\Delta }$ 在上与LGL配置点对应的点。

3. LST-LGNI格式

引进通量 $u=\nabla p$ 则

$\nabla \times \left(u\right)|{\text{Γ}}_N=0,{\gamma }_{\tau }u|{\text{Γ}}_D=0,$

其中 ${\gamma }_{\tau }u=n\times u$ ，n和 ${\gamma }_{\tau }$ 分别对应边界 ${\Gamma }_D$ 的单位法向量和单位切向量。如文 [1] 得到一阶方程组为

$\{\begin{array}{l}u-\nabla p=0,\text{in}{\Omega }_{\Delta }\\ -\nabla \cdot u+b\cdot u+cp=f\left(x,y\right),\text{in}{\Omega }_{\Delta }\\ \nabla \times u=0,\text{in}{\Omega }_{\Delta }\end{array}$ (3.1)

其中边界条件为

$p|{\text{Γ}}_D=0,n\cdot u|{\text{Γ}}_N=0,{\gamma }_{\tau }u|{\text{Γ}}_D=0.$

设

$V:=\left\{v\in H^1\left(\Omega \Delta \right):v|{\Gamma }_D=0\right\},$

$W:=\left\{v\in {\left[H^1\left(\Omega \Delta \right)\right]}^2:n\cdot v|{\Gamma }_D=0,\tau v|{\Gamma }_D=0\right\}.$

定义系统(3.1)在 $L^2$ 范数下的最小二乘函数为：

$F\left(v,p;f\right)={\Vert f+\nabla \cdot v-b\cdot v-cq\Vert }_{{\Omega }_{\Delta }}^2+{\Vert v-\nabla q\Vert }_{{\Omega }_{\Delta }}^2+{\Vert \nabla \times v\Vert }_{{\Omega }_{\Delta }}^2,$ (3.2)

其中 $\left(v,q\right)\in W\times V$ 。即，对应的最小值问题为：求 $\left(u,p\right)\in W\times V$ 满足

$F\left(u,p;f\right)=\underset{\left(v,q\right)\in W\times V}{\inf }F\left(v,q;f\right).$

因此，上述对应的变分格式为：求 $\left(u,p\right)\in W\times V$ 满足

$a\left(u,p;v,q\right)=f\left(v,q\right),\forall \left(v,q\right)\in W\times V,$ (3.3)

其中

$\begin{array}{l}a\left(u,p;v,q\right)=\left(\nabla \cdot u-b\cdot u-cp,\nabla \cdot v-b\cdot v-cp\right){\Omega }_{\Delta }\\ +\left(u-{\nabla }_p,v-{\nabla }_q\right)+\left(\nabla \times u,\nabla \times v\right){\Omega }_{\Delta },\\ f\left(v,q\right)=-\left(f,\nabla \cdot v-b\cdot v-cq\right){\Omega }_{\Delta }.\end{array}$

设 ${\mathbb{Q}}_N\left({\Omega }_{\Delta }\right)$ 为三角形区域 ${\Omega }_{\Delta }$ 上对应有限维空间 ${\hat{\mathbb{Q}}}_N$ 的所有次数不超过N多项式的空间。设

$W_N={\left[{\mathbb{Q}}_N\left({\Omega }_{\Delta }\right)\right]}^2\cap W,V_N={\mathbb{Q}}_N\left({\Omega }_{\Delta }\right)\cap V.$ (3.4)

最小二乘函数(3.2)离散形式为：

$F_N\left(v,q;f\right)={\Vert I_{N}f+\nabla \cdot v-I_N\left(b\cdot v+cq\right)\Vert }_{{\Omega }_{\Delta },N}^2+{\Vert I_{N}v-\nabla q\Vert }_{{\Omega }_{\Delta },N}^2+{\Vert \nabla \times I_{N}v\Vert }_{{\Omega }_{\Delta },N}^2,$

其中 $\left(v,q\right)\in W_N\times V_N$ ， $I_N$ 为Legendre-Gauss-Lobatto (LGL)插值算子，或者Chebyshev-Gauss-Lobatto (CGL)配插值算子。

类似地，得到离散的最小问题为：求 $\left(u_N,p_N\right)\in W_N\times V_N$ 使得

$F_N\left(u_N,p_N;f\right)=\underset{\left(v,q\right)\in W_N\times V_N}{\inf }{F}_N\left(v,q;f\right).$ (3.5)

因此，最小二乘三角形Legendre Galerkin数值积(LST-LGNI)格式为：求 $\left(u_N,p_N\right)\in W_N\times V_N$ 使得

$a_N\left(u_N,p_N;v,q\right)=f_N\left(v,q\right),\forall \left(v,q\right)\in W_N\times V_N,$ (3.6)

其中

$\begin{array}{l}{a}_N\left(u_N,p_N;v,q\right)=\left(\nabla \cdot u_N-I_N\left(b\cdot u_N+c{p}_N\right),\nabla \cdot v-I_N\left(b\cdot v+cq\right)\right){\Omega }_{\Delta },N\\ +\left(u_N-\nabla p_N,v-\nabla q\right){\Omega }_{\Delta },N+\left(\nabla \times u_N,\nabla \times v\right){\Omega }_{\Delta },N,\end{array}$

$f_N\left(v,q\right)=-\left(I_{N}f,\nabla \cdot v-I_N\left(b\cdot v+cq\right)\right){\Omega }_{\Delta },N.$

4. 数值算例

为了验证格式(3.6)的谱精度和有效性，本节给出一些数值算例且与一些三角单元谱方法进行比较。下面定义离散的 $L^2$ 误差为：

$E_N={\left(\underset{l,m=0}{\overset{N}{{\displaystyle \sum }}}{\left(u_N\left(x_l,y_m\right)-u\left(x_l,y_m\right)\right)}^2{w}_{lm}\right)}^{1/2},$

其中 ${\left\{\left(x_l,y_m\right)\right\}}_{l,m=0}^N$ 和 ${\left\{w_{lm}\right\}}_{l,m=0}^N$ 分别为LGL节点和其对应的数值积分权函数。

例4.1：考虑问题(4.1)， $b=\left(0,0\right)$ ，以及 $c=0$ 。精确解为

$p\left(x,y\right)=xy\left(e^{\left(x+y\right)}-e\right),$ (4.1)

$p\left(x,y\right)=\frac{xy\left(1-x-y\right)}{\left(x+1/10\right)\left(y+1/10\right)}.$ (4.2)

文 [17] 中采用Chebyshev配置方法对该算例进行计算，这里则用LST-LGNI (3.6)计算。表1中分别给出了文 [18] 和我们的格式(3.6)的关于L²-范数的数值误差，其中误差的定义为 ${\bar{E}}_N=\frac{E_N}{N}$ 以及 $w_{lm}=1$ 。数值结果表明，该方法具有高阶谱精度且文 [18] 中给出的数值结果类差不多。

Table 1. Error in L2-norm for example 4.1

表1. 例4.1中的L²-范数误差

例4.2：考虑问题(4.4)，其中 $b=\left(0,0\right),c=1$ 。精确解分别为(4.1)和图1中给出了最大值点误差和离散的L²-误差且数值结果显示误差具有高阶谱精度。表2中给出数值表明该方法的谱精度优于 [19] 方法。

Figure 1. Maximum point error and L²-norm error; On the left hand side: Error of the exact solution (4.1), on the right hand side: Error of the exact solution (4.3), where $\alpha =3.2$ , $\beta =2$

图1. 最大值点误差和L²-范数误差。左侧：精确解为(4.1)的误差情况；右侧：精确解(4.3)的误差情况，其中 $\alpha =3.2$ 和 $\beta =2$

例4.3：考虑问题 [11]

$\{\begin{array}{l}-\nabla \cdot \left(\left(x+2\right)\nabla p\right)+b\cdot \nabla p+\left(x+y\right)p=f\left(x,y\right),\left(x,y\right)\in {\Omega }_{\Delta },\\ p=0,\left(x,y\right)\in {\Gamma }_D,\\ n\cdot \nabla p=0,\left(x,y\right)\in {\Gamma }_N.\end{array}$ (4.4)

精确解如文 [11] 中给出为

Table 2. Error in L2-norm for example 4.2

表2. 例4.2中的L²-范数误差

$p\left(x,y\right)=e^{x+y-1}\sin \left(3y\left(y-\sqrt{3}\frac{x}{2}+\frac{\sqrt{3}}{4}\right)\right).$ (4.5)

为了验证我们的方法LST-LGNI (3.6)对变系数为多项式情况也是有效的。在该算例中将我们的数值结果与文 [11] 的三角形Legendre谱方法的进行比较。表3中分别给出文 [11] 中三角谱方法与LST-LGNI (3.6)的数值结果。由此注意到格式(3.6)谱精度比 [11] 中给出的要稍好一些。

Table 3. Error in L2-norm for example 4.3

表3. 例4.3中的L²-范数误差

5. 总结

本文研究了最小二乘Legendre Galerkin的数值积分的方法。该方法结合最小二乘，离散时采用Legendre Galerkin数值积分处理，使得格式计算方便且导出的代数方程系数具有对称正定的特点，求解时允许采用如共轭梯度法。数值实验验证该方法具有高阶谱精度。更有意义是考虑复杂区域上的微分方程的求解，类似文 [14] 进一步将结合三单元与矩形元发展高效的最小二乘谱元格式。

基金项目

国家自然科学基金资助项目(No.11701119)，大学生创新训练资助项目(No.201710595033)，广西自然科学基金资助项目(No.2017GXNSFBA198053)，广西混杂计算与集成电路设计分析重点实验室开放课题资助(No.HCIC201607)。

参考文献

NOTES

*通讯作者