自从Reed和Hill在文献[1]中处理双曲方程时首次提出间断Galerkin方法(简称DG方法)以来, DG方法作为求解双曲问题以及类双曲问题的有效方法而蓬勃发展.随后, Babuska等在文献[2]中提出了求解椭圆和抛物问题的间断Galerkin有限元方法.还有许多的DG方法被研究工作者提出, 参见文献[3-12]等.这些方法一般被称作内部惩罚方法.此外, Cockburn和Shu^[13]提出了局部间断的Galerkin有限元方法(LDG方法); Oden和Bauman^[14-15]在处理某些扩散问题时又提出了其他几种形式的DG方法.

然而, 上面提及的方法在应用时需事先确定稳定化参数, 而这些参数的确定是极其不易的.基于文献[16-17]的思想, 本文提出一种稳定化的$hp$间断Galerkin方法, 给出满足局部守恒性质的新变分格式, 且该方法的稳定化参数易确定.

本文其余部分结构如下, 第二节引入模型问题并给出稳定化间断Galerkin方法.第三节推导出最优阶的误差估计, 并分析该数值方法的收敛性.第四节通过数值算例验证理论结果的正确性.最后, 给出结论总结本文的主要结果.

考虑如下反应扩散问题

$ - \nabla \cdot (K(x)\nabla u) + u = f\textrm{in}\ \Omega, $

(2.1)

$ u = 0\textrm{on}\ \partial \Omega, $

(2.2)

其中$\Omega \subset {\Bbb R}^d(d=2, 3)$为有界开区域, 其边界$\partial \Omega$是Lipschitz连续的. $f\in L^d(\Omega)$, 以及$K(x)\in L^\infty(\Omega)$使得$0 < K_1 \leq K(x) \leq K_2$.

需要说明的是, 问题(2.1)-(2.2)解的正则性可能很低, 因此, 对于其变分问题的弱解, 需要定义合适的解空间.为此, 我们先引入一些记号.记$\{P_h\}$是对$\Omega$的正则剖分族, $P_h$中的元素$E$均为无交的开集, 即$\Omega = \textrm{int} \big(\bigcup\limits_{E \in P_h} \bar{E}\big)$.记$h_E=\textrm{diam}(E)$, $h=\max\limits_{E \in P_h}h_E$. $P_h$中所有边界构成的集合记为$\varepsilon_h = \{\gamma_k\}, k=1, \ldots, N_{edge}$, 其中$N_{edge}$为剖分$P_h$中边界的总数.而剖分$P_h$的所有内部单元边界的集合$\Gamma_{\textrm{int}}$, 记为$\Gamma_{\textrm{int}} = \bigcup\limits_{k=1}^{N_{edge}}\gamma_k \backslash \partial \Omega$.

定义单元边界$\gamma_e \in \Gamma_\textrm{int}$上的跳跃(Jump)和平均(Average)分别为

$ [v]=v|_{\gamma_e \subset \partial E_i}-v|_{\gamma_e \subset \partial E_j}, \;\;\langle v\rangle =\frac12(v|_{\gamma_e \subset \partial E_i} + v|_{\gamma_e \subset \partial E_j}), \ i>j, $

(2.3)

其中$\gamma_e=\textrm{int}(\partial E_i \cap \partial E_j)$是二维相邻单元的公共边(或三维相邻单元的公共面), 如图 1和图 2所示.

为给出问题(2.1)–(2.2)的变分格式, 进而求得弱解, 首先引入分片Sobolev空间

$ {\cal M}(P_h)=\big\{v \in L^2(\Omega): v|_E \in H(\Delta, E), [\nabla v \cdot { \boldsymbol{n}}] \in L^2(\Gamma_{\textrm{int}}), \forall E \in P_h\big\}, $

(2.4)

其中$H(\Delta, E)=\{v \in L^2(E): \nabla \cdot \nabla v \in L^2(E)\} \subset H^1(E)$, $M(P_h)$的范数定义为

$ |||v|||^2 = \sum\limits_{E \in P_h}\|v\|^2_{H^1(E)} +\sigma\frac{h^\lambda}{p^\zeta}\| [K(x)\nabla v \cdot {\boldsymbol{n}}]\|^2_{L^{2}(\Gamma_{\textrm{int}})} {}, $

(2.5)

这里$p$是多项式空间中近似多项式次数的最小值, 参数$\lambda$和$\zeta$均为正常数.

记$V$是${\cal M}(P_h)$的关于范数$|||\cdot|||$的完备化空间.基于前面的准备工作, 现给出问题(2.1)–(2.2)的变分问题:求$u \in V$使得

$ B(u, v)=L(v), \;\;\forall v \in V, $

(2.6)

其中$B(u, v)$和$L(v)$分别为

$ B(u, v) = \sum\limits_{E \in P_h}\Big\{\int_{E}\big(K(x)\nabla u \cdot \nabla v + uv\big){\rm d}x {}\nonumber\\ \;\;\;\;\;\;\;\;\;\;\;\;- \int_{\partial E}\big(v(K(x)\nabla u \cdot {\boldsymbol{\mu}}) -(K(x)\nabla v \cdot {\boldsymbol \mu})u\big){\rm d}s\Big\} {} \nonumber\\ \;\;\;\;\;\;\;\;\;\;\;\;+ \int_{\Gamma_{\textrm{int}}}\big(\langle v\rangle[K(x)\nabla u \cdot {\boldsymbol{n}}] -\langle u\rangle[K(x)\nabla v \cdot {\boldsymbol{n}}]\big){\rm d}s {} \nonumber\\ \;\;\;\;\;\;\;\;\;\; \;\;+ \int_{\Gamma_{\textrm{int}}} \sigma \frac{h^{\lambda}}{p^{\zeta}} [K(x)\nabla u \cdot {\boldsymbol{n}}][K(x)\nabla v \cdot {\boldsymbol n}]{\rm d}s, {} $

(2.7)

$ L(v) = \int_{\Omega}fv {\rm d}x. $

(2.8)

注2.1  与文献[17-18]中方法的不同之处在于, 本文的方法可以处理一类解具有较低正则性的反应扩散方程.本文定义的双线性形式$B(u, v)$也不同于文献[19]中给出的双线性形式, 本文中考虑的跳跃项为$[\nabla v\cdot { \boldsymbol{n}}]$而文献[19]中的跳跃项为$[v]$.

设$\{F_E\}$是定义在剖分$P_h$上的一族可逆仿射变换: $F_E: \hat{E} \rightarrow E, x=F_E(\hat{x})$, 其中$\hat{E}$为参考单元, $\hat{x}$为参考单元上的点.

定义如下有限元空间

$ V_h=\{v \in L^2(\Omega)|\, v|_E=\hat{v}\circ F^{-1}_E, \hat{v} \in P_{p_E}(\hat{E}), \forall E \in P_h\}. $

(2.9)

容易验证$V_h\subset V$, 其中$p_E$为$E$上分片多项式的次数.

基于(2.6)式, 我们给出离散变分问题:求$u_h \in V_h$使得

$ B(u_h, v_h) = L(v_h), \forall v_h \in V_h, $

(2.10)

其中双线性形式$B(\cdot, \cdot)$由(2.7)式给出.

对于双线性形式$B(\cdot, \cdot)$, 我们有如下结论.

定理2.1  存在正常数$M>0$, 使得

$ |B(u_h, v_h)|\leq M|||u_h|||\, |||v_h|||, \ \forall u_h, v_h \in V_h. $

(2.11)

证  详细证明, 参见文献[16].

定理2.2  存在常数$\gamma_h = \gamma(\sigma, h, p)>0$使得

$ \sup\limits_{v_h \in V_h\backslash\{0\}} \frac{|B(u_h, v_h)|}{|||v_h|||} \geq \gamma_h |||u_h|||, \;\;\;\;\forall u_h \in V_h. $

(2.12)

证  由上确界的定义, 可得

$ \sup\limits_{v_h \in V_h\backslash\{0\}} \frac{|B(u_h, v_h)|}{|||v_h|||} \geq \frac{B(u_h, u_h)}{|||u_h|||}. $

(2.13)

应用(2.7)式并取$C=\min(K_1, 1)$, 有

$ B(u_h, u_h) =\sum\limits_{E \in P_h}\int_E\Big(K(x)|u_h|^2+u_h^2\Big) {} +\sigma\frac{h^\lambda}{p^\zeta}\| [K(x)\nabla u_h \cdot {\boldsymbol{n}}]\|^2_{L^{2}(\Gamma_{\textrm{int}})} {}\nonumber\\ \geq C\sum\limits_{E \in P_h}\|u_h\|_{H^1(E)}^2+K_1\sigma\frac{h^\lambda}{p^\zeta}\| [\nabla u_h \cdot { \boldsymbol{n}}]\|^2_{L^{2}(\Gamma_{\textrm{int}})} {}\nonumber\\ \geq C|||u_h|||. $

(2.14)

联合(2.13)和(2.14)式, 即证得定理成立.

注2.2  定理2.1和定理2.2表明双线性形式$B(\cdot, \cdot)$满足连续性和强制性条件, 因此根据有限元理论^{[5, 20]}, 变分问题(2.10)解存在唯一.

为了给出该方法的最优误差估计, 首先引入一系列插值, 令$\tilde{\pi}^{E}_{hp}: H^{\gamma_{k}}(E) \rightarrow P_{p_{E}}(E)$满足

$ \int_{\gamma \subset \partial E}\nabla (\varphi - \tilde{\pi}^{E}_{hp}(\varphi)) \cdot {\boldsymbol{\mu}}{\rm d}s = 0. $

(3.1)

对单元$E$将局部插值算子$\tilde{\pi}^{E}_{hp}$零延拓到整个剖分$P_h$上, 定义全局插值$\tilde{\Pi}_{hp}:V \rightarrow V_h$, 即

$ \tilde{\Pi}_{hp}(u)=\sum\limits_{E \in P_h} \tilde{\pi}^{E}_{hp}(u|_E), u \in V. $

(3.2)

对于上述的插值算子, 有如下的不等式成立(参见文献[12])

$ \|\varphi - \tilde{\pi}^{E}_{hp}(\varphi)\|_{L^2(E)} \leq C\frac{h^{\mu_k}}{p^{\gamma_k-3/2}_{E}}\|\varphi\|_{H^{\gamma_k}(E)}, $

(3.3)

$ \|\nabla (\varphi - \tilde{\pi}^{E}_{hp}(\varphi))\|_{L^2(E)} \leq C\frac{h^{\mu_k-1}}{p^{\gamma_k-3/2}_{E}}\|\varphi\|_{H^{\gamma_k}(E)}, $

(3.4)

其中$\mu_k = \min\{p_E+1, \gamma_k\}, \ p_E \geq 1, \ \gamma_k\geq 2$.

此外, 迹定理^[20]成立

$ \|\nabla w \cdot {\boldsymbol{\mu}}\|^2_{L^2(\partial E)} \leq C \Big\{ \frac{1}{h_E}\|\nabla w \|^2_{L^2(E)}+\|\nabla w \|_{L^2(E)} \|\nabla^2 w \|_{L^2(E)}\Big\}, $

(3.5)

其中$C>0$是不依赖与$h_E$的常数.

由有限维空间范数的等价性, 我们有

引理3.1  存在依赖$\sigma$的常数$C$, 使得

$ |||v_h-\bar{v}_h|||\leq C(\sigma) p\|v_h\|_{H^1(P_h)}, \ \forall v_h \in V_h, $

(3.6)

其中$\bar{v}_h = \sum\limits_{E \in P_h} \bar{v}_h|_E, \bar{v}_h|_E = \frac{1}{|E|}\int_E v_h {\rm d}x$.

记$e_h = u-u_h = u-\tilde{\Pi}_{hp}(u) + \tilde{\Pi}_{hp}(u)-u_h=\eta_h + \xi_h$, 这里$\eta_h=u-\tilde{\Pi}_{hp}(u), \ \xi_h=\tilde{\Pi}_{hp}(u)-u_h$, 那么

$ B(e_h, v)=L(v)-B(u_h, v)=R^h(v), \forall v \in V, $

(3.7)

其中$R^h:V \rightarrow {\Bbb R}$称为残差泛函.由于残差在$V_h$上的正交性, 有

$ B(e_h, v_h)=0, \forall v_h \in V_h. $

(3.8)

定理3.1  设$u \in H^2(\Omega) \cap V$, 则存在不依赖于$h$和$p$的常数$C>0$, 使得

$ |||u-\tilde{\Pi}_{hp}(u)|||\leq C\frac{h^{\mu^{\ast}}}{p^{\gamma^{\ast}}} \|u\|_{H^{\gamma_k}(P_h)},\;\;\; \gamma_k \geq 2, \ p \geq 1, \ $

(3.9)

其中$\mu^{\ast}=\min \big\{\mu-1, \mu-\frac32+\frac\lambda 2 \big\}$, $\gamma^{\ast}=\min \big\{\gamma-\frac32, \gamma-\frac32+\frac{\zeta}2\big\}$, $\mu=\min \big\{p+1, \mu_k\big\}$和$\gamma=\min\limits_{E \in P_h}(\gamma_k)$.

证  由(3.5)式, 得

$ |||\eta_h|||^2 \leq \sum\limits_{E \in P_h} \|\eta_h\|^2_{H^1(E)} +K_2^2\sigma \frac{ h^{\lambda}}{p^\zeta}\|[\nabla \eta_h\cdot {\boldsymbol{\mu}}]\|^2_{L^{2}(\Gamma_{\textrm{int}})}\nonumber\\ \leq \sum\limits_{E \in P_h} \|\eta_h\|^2_{H^1(E)} +C K_2^2\sigma \frac{ h^{\lambda-1}}{p^\zeta}\|[\nabla \eta_h\cdot {\boldsymbol{\mu}}]\|^2_{L^{2}(\Gamma_{\textrm{int}})}. $

(3.10)

应用(3.10)式和(3.3)–(3.4)式, 即得(3.9)成立.定理得证.

定理3.2  设$u \in H^2(\Omega) \cap V$是问题(2.6)的解, 且$u_h \in V_h$满足(2.10)式.若稳定化参数$\lambda \geq 1$, $\zeta \geq 0$, 则存在常数$C\geq 0$, 使得

$ \|u-u_h\|_{H^1(P_h)} \leq C \frac{h^{\mu-1}}{p^{\gamma-5/2}} \|u\|_{H^{\gamma_k}(P_h)}, p \geq 1, $

其中$\mu = \min \{p+1, \gamma\}$.

证  应用三角不等式, 有

$ \|u-u_h\|_{H^1(P_h)} \leq \|\eta_h\|_{H^1(P_h)}+\|\xi_h\|_{H^1(P_h)}. $

(3.11)

由(2.11), (2.12), (3.8)式以及引理3.1, 得

$ \|\xi_h\|^2_{H^1(P_h)} \leq B(\eta_h, \xi_h)=B(\eta_h, \xi_h-\bar{\xi_h}) +B(\eta_h, \bar{\xi_h})\nonumber\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\leq C|||\xi_h-\bar{\xi_h}||| |||\eta_h|||+B(\eta_h, \bar{\xi_h})\nonumber\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\leq C(\sigma)p\|\xi_h\|_{H^1(P_h)} |||\eta_h|||+B(\eta_h, \bar{\xi_h}), $

(3.12)

其中$\bar{\xi_h}$是$\xi_h$的分片平均.

由(2.7)式和柯西-施瓦兹不等式, 有

$ B(\eta_h, \bar{\xi_h})=\sum\limits_{E \in P_h}\Big(\int_E \eta \bar{\xi_h}{\rm d}x -\int_{\partial E}\bar{\xi_h}K(x)\nabla \eta_h \cdot {\boldsymbol{\mu}}{\rm d}s\Big)\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+\int_{\Gamma_\textrm{int}} \langle \bar{\xi_h} \rangle [K(x)\nabla \eta_h \cdot { \boldsymbol{n}}]{\rm d}s\nonumber\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;=\sum\limits_{E \in P_h} \int_E \eta_h \bar{\xi_h}{\rm d}x\nonumber\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\leq \|\eta_h\|_{L^2(\Omega)}\|\bar{\xi_h}\|_{L^2(\Omega)}\nonumber\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\leq C\|\eta_h\|_{L^2(\Omega)}\|\xi_h\|_{H^1(P_h)}. $

(3.13)

联合(3.11), (3.13)和(3.12)式, 可得

$ \|u-u_h\|_{H^1(P_h)} \leq C p |||\eta_h|||. $

(3.14)

应用定理3.1和(3.14)式, 得

$ \|u-u_h\|_{H^1(P_h)} \leq C \frac{h^{\mu^{\ast}}}{p^{\gamma^{\ast}-1}} \sqrt{\sum\limits_{E \in P_h} \|u\|^2_{H^{\gamma_k}(E)}}, $

(3.15)

其中$\mu^{\ast}=\mu -1$和$\gamma^{\ast} = \gamma -\frac{3}{2}$.定理证毕.

定理3.3  设$u \in H^2(\Omega) \cap V$是问题(2.6)的真解, 数值解$u_h \in V_h$满足(2.10)式.若稳定化参数$\lambda \geq 1$和$\zeta \geq 0$, 则存在不依赖于$h$和$p$的常数$C>0$, 使得

$ \|u-u_h\|_{L^2(\Omega)} \leq C \frac{h^{\mu}}{p^{\gamma-5/2}} \|u\|_{H^{\gamma_k}(P_h)}, p \geq 1. $

(3.16)

证  为证明$L^2$范数下的误差估计, 应用Aubin-Nitsche技巧(参见文献[21]).为此, 引入如下问题

$ - \nabla \cdot (K(x)\nabla \phi) + \phi = e_h\textrm{in}\ \Omega, $

(3.17)

$ \phi = 0\quad\ \textrm{on}\ \partial \Omega. $

(3.18)

对于上述问题(3.17)-(3.18), 有正则性结果

$ \|\phi\|_{H^2(\Omega)} \leq C\|e_h\|_{L^2(\Omega)}. $

(3.19)

由(3.17)式知

$ \|e_h\|^2_{L^2(\Omega)}=\int_{\Omega}(- \nabla \cdot (K(x)\nabla \phi) + \phi)e_h{\rm d}x\nonumber\\ =\sum\limits_{E \in P_h}\int_E (K(x)\nabla\phi \cdot \nabla e_h+\phi e_h){\rm d}x-\int_{\Gamma_{\rm int}}\langle K(x)\nabla \phi \cdot { \boldsymbol{\mu}}\rangle [e_h]{\rm d}s. $

(3.20)

联合(3.8)和(3.20)式, 得

$ \|e_h\|^2_{L^2(\Omega)}=\sum\limits_{E \in P_h}\int_E (K(x)\nabla(\phi-v_h) \cdot \nabla e_h+ (\phi-v_h)e_h){\rm d}x{}\nonumber\\ {} -\int_{\Gamma_{\rm int}}\langle K(x)\nabla(\phi-v_h) \cdot {\boldsymbol{\mu}} \rangle [e_h]{\rm d}s -\int_{\Gamma_{\rm int}}\langle K(x)\nabla e_h \cdot {\boldsymbol{n}}\rangle [v_h]{\rm d}s{}\nonumber\\ {} -\int_{\Gamma_{\rm int}} \sigma\frac{h^\lambda}{p^\zeta} [K(x)\nabla e_h\cdot {\boldsymbol{n}}][K(x)\nabla v_h\cdot {\boldsymbol{n}}]{\rm d}s{}\\ = A_1 + A_2, $

(3.21)

其中$A_1=\sum\limits_{E \in P_h}\int_E(K(x)\nabla(\phi-v_h)\cdot \nabla e_h + (\phi-v_h)e_h){\rm d}x$, $A_2$为(3.21)式中除去$A_1$余下的项.

取$v_h = \tilde{\phi}$为$\phi$的次数为$p_E$的连续插值, 应用柯西-施瓦兹不等式, 有

$ |A_1| \leq C\|\phi\|_{H^1(\Omega)} |||e_h||| \leq Ch\|\phi\|_{H^2(\Omega)} |||e_h|||, $

(3.22)

$ |A_2| \leq |B(e_h, \phi-\tilde{\phi})| \leq M |||e_h|||\, |||\phi-\tilde{\phi}|||. $

(3.23)

利用(3.19), (3.22), (3.23)式以及定理3.2, 可得

$ \|e_h\|^2_{L^2(\Omega)} \leq C \frac{h^{\mu}}{p^{\gamma-5/2}} \|e_h\|_{L^2(\Omega)} \|u\|_{H^{\gamma_k}(P_h)}, $

(3.24)

即得(3.3)式成立.定理得证.

本文仅就二维问题给出数值实验.在单位正方形$\Omega = (0, 1)\times(0, 1)$上考虑问题(2.1)-(2.2), 取$K(x, y) = xy$, 真解取为$u(x, y) = xy(1-x)(1-y)$.稳定化参数取为$\lambda=\zeta=1$.记$N=-\ln(h)/\ln(2)$, $E_N$为相应范数的误差, $e_h=E_N/E_{N-1}$, $r=\ln e_h/\ln2$, 下列表格给出不同范数、不同网格密度的误差结果.

表 1-表 4分别给出了应用$P_1$元和$P_2$元在$L^2$范数和$H^1$范数下的误差和收敛率.对于应用高阶元-$P_3$元和$P_4$元的误差结果, 本文不再列出.

下面图 3和图 4分别给出了$L^2$范数和$H^1$范数下的收敛率变化情况.

作为比较, 图 5和图 6分别给出了真解和数值解的等值线图.

本文提出了求解反应扩散问题的一种新的绝对稳定$hp$间断Galerkin有限元方法.由注2.1知本文的方法是新的, 此外该文的创新之处还有:稳定化参数$\lambda \geq 1$, $\zeta \geq 0$时, 得到了最优的收敛阶; 定义了简单而自然的与网格有关的范数; 对于$p\geq 1$, 定理3.2和定理3.3都成立.而文献[17]中, $p$须满足$p\geq 2$.最后, 数值算例验证了理论结果的正确性.