扩散反应方程的一种协调间断有限元法

顾子兵 胡朝浪 杨荣奎 冯民富

摘 要    ：对于定常扩散反应方程的原始混合变分形式，本文基于间断有限元法和最小二乘法思想给出了一种新的协调间断有限元格式，并将其推广应用于非定常扩散反应方程，进而建立了全离散协调间断有限元格式. 该格式不仅能避免LBB条件，而且适用于多边形网格. 在能量范数下，本文得到了原始变量和通量的最优误差估计.最后，针对定常和非定常的扩散反应方程，本文分别以一般情形和对流占优情形下的数值算例验证了方法的有效性.

关键词 ：协调间断有限元； 原始混合形式； 最小二乘法

中图分类号 ： O241.82 文献标识码 ：A DOI ：  10.19907/j.0490-6756.2023.041004

A conforming discontinuous finite element  method for diffusion-reaction equations

GU Zi-Bing， HU Chao-Lang， YANG Rong-Kui， FENG Min-Fu

（School of Mathematics， Sichuan University， Chengdu 610064， China）

For the primal-mixed form of the steady diffusion-reaction equation，  a new conforming discontinuous finite element scheme is proposed by using the idea of discontinuous finite element method and the least-square method. Then  the method is extended to the case of unsteady diffusion-reaction equation and  a fully discrete conforming discontinuous finite element scheme is given. This scheme avoids the LBB condition and can be applied to polygon mesh. Under the energy norm， the optimal convergence order error estimates for the primary and flux variables are established. Finally， for the steady and unsteady equations， some numerical examples in general and convection dominated cases are given to verify the effectiveness of the method.

Conforming discontinuous element; Primal-mixed form; Least-square method

（2010MSC 65M60）

1 引 言

本文考虑如下扩散反应问题：求函数 u=u（x） ，使得

-  （1）

其中，Ω是  Euclid Math TwoRA@

d d=2，3  中的一个多边形区域，其Lipschitz边界为   Ω， A 为正定矩阵， c≥0 . 该问题常被用于描述温度、压力及化学物质的浓度等物理量的演化.

有限元法是求解问题（1）的一种有效手段. 用经典的有限元法可以直接求得问题（1）的离散解. 在实际问题中， u 的梯度

u 往往有重要的物理意义，如在流动问题中

u 表示流体的速度. 为求 u 和

u， 2013年Wang和Ye  ［1］对二阶椭圆问题提出了弱有限元法. 该方法区别于间断有限元法的地方在于它将一个不连续的函数 v 在单元上的取值分别为内部  v  0 和边界  v  b ，在构造有限元时分别在内部和边界上构造函数，用弱梯度

w· 代替梯度

. 弱有限元法具有和间断有限元法相似的优点，且形函数选取及网格剖分多样. 2014年，Wang和Ye  ［2］针对带Dirichlet边界的扩散方程的对偶混合形式提出了弱Galerkin混合有限元方法.该方法对主要变量和通量都有非常精确地近似，且可以应用于多边形和多面体网格. 2017年，Mu等  ［3］提出了最小二乘的弱有限元方法，该方法可以被用于奇异摄动扩散反应方程和各向异性、不均匀性多孔介质中流体的流动问题，具有较好的鲁棒性和灵活性. 值得注意的是，上述工作都是基于方程的对偶混合变分形式提出的方法.

2014年，Wang等  ［4］提出了修正的弱有限元法.该方法将弱有限元方法中的边界函数  v  b 替换为  v  ，通过单元内部函数 v 确定边界函数  v  ， 减少了边界自由度，提高了计算效率. 2020年，Ye等  ［5-8］提出了协调间断元方法.该方法属于弱有限元方法，边界函数仍然为  v  . 此外，该方法还具备协调有限元和间断有限元方法的特征，形式简洁，单元形函数采用间断函数，且可以通过调整有限元空间将该方法应用于多边形和多面体网格.

基于上述工作，本文考虑方程（1）的原始混合变分形式. 引入通量 σ=-A

u 后，方程（1）的求解等价于如下问题：

求  u，σ ∈ H  1  0（ Ω ）×   L  2（ Ω ）   d ，使得

A  -1σ，τ +

u，τ =0， τ∈   L  2（ Ω ）   d， （σ，-

v）+（cu，v）=〈f，v〉， u∈ H  1  0（ Ω ）   （2）

对问题（2），本文采用协调间断有限元方法结合最小二乘法  ［3，9］，给出了一种新的协调间断有限元格式. 通过引入最小二乘的残差项，该方法不但避免了LBB条件  ［10］的限制， 而且可被应用于多边形网格.

本文结构如下. 第2节针对定常扩散反应方程给出了一种协调间断有限元方法的数值格式，并给出理论分析. 在此基础上，本文也给出一般情况和反应占优时的数值算例，验证了该格式的有效性. 在第3节中，本文将第2节中的协调间断有限元格式推广应用于非定常扩散反应方程，给出了一种全离散协调间断有限元格式，并通过一般情况和反应占优时的数值算例验证了该格式的有效性. 第4节总结本文所得结果.

2 定常扩散反应方程情形

2.1 协调间断有限元格式

我们以  W  m，p（ Ω ） ，  W  m，p  0（ Ω ） 表示Ω上的 m -阶Sobolev空间. 特别地，当 p =2时，

H  m（ Ω ）=W  m，2（ Ω ），     H  m  0（ Ω ）=W  m，2  0（ Ω ），

·，·   m ，   ·   m ，   ·   m 分别为  H  m（ Ω ） 的内积、范数和半范数.

令  T  h 为区域  Ω  的有限元剖分，该剖分由凸且正则的多边形或多面体组成. 记  F  h 为  T  h 的所有边或面的集合，  F  0  h= F  h＼  Ω  表示所有内边的集合. 为了简便，分别定义  T  h 和   T  h 上的内积为

v，w     T   h=∑  T∈ T   h  （v，w）   T=∑  T∈ T   h ∫   Tvw d x，

v，w      T   h=∑  T∈ T   h   v，w     T=∑  T∈ T   h ∫    Tvw d s.

为了简便表示函数的弱梯度，用 v 而不用  v  h 表示有限元空间  V  h 的任意函数.

对每个单元 T∈ T  h ，记  h  T 为其直径，  T  h 的网格大小为 h=  max   T∈ T  h    h  T . 定义弱函数 v ，

v= v， 在单元T内部，

{v}， 在边界 T上.

设整数 k≥1 ，记  P  k（T） 为定义在 T 上次数不超过 k 的多项式空间， ［ P  k（T） ］  d 为有限维向量多项式空间.  v 的有限元空间  V  h 可如下定义：

V  h={v∈ L  2（ Ω ）：v |  T∈［ P  k（T）］， T∈ T  h，

v |    Ω =0}.

通量 σ 的有限元空间  W  h 可定义如下：

W  h={τ∈ ［ L  2（ Ω ）］  d：τ |  T∈ ［ P  k-1（T）］  d，

T∈ T  h}.

接下来，定义 v 的弱梯度

wv . 令  T  1 和  T  2 为  T  h 中包含公共边 e∈ F  h 的两个单元. 对于 e∈ F  h 和  v∈ V  h 的跳量  v  定义为

v =v，e?  Ω ， v =v |   T  1-v |   T  2，e∈ F  0  h.

这里，  T  1 和  T  2 的顺序不重要. 对于 e∈ F  h 和 v∈ V  h ， v 的均值  v  定义为

v =v，e?  Ω ，{v}= 1 2 （v |   T  1+v |   T  2），

e∈ F  0  h.

对于函数 v∈ V  h ，它的弱梯度

wv 是一个分片向量多项式，即

wv∈∏  T∈ T  h［ P  j（T） ］  d，j>k . 在每个单元 T 上，对任意 q∈［ P  j T  ］  d ，

wv 满足以下方程：

（

wv，q ）  T=（v，-

·q ）  T+ 〈 v ，q·n〉   T   （3）

它的弱散度

w·v 是一个分片多项式，即

w·v∈∏   T∈ T   h P   k-1（T） . 对任意 的τ∈［ P  k-1（T） ］  d ，在每个单元 T 上

w·v 满足以下方程：

（

w·v，τ ）  T=（v，-

τ ）  T+ 〈{v}·n，τ〉   T.

注  在方程（1）～（3）中， j 的选取依赖于多边形的边数. 对于三角形网格， j 可以取 k+1 .通常情况下， j=n+k-1 ，其中 n 为多边形的边数  ［11］.

问题（2）等价于如下弱形式.

求 （u，σ）∈ H  1  0（ Ω ）×   L  2（ Ω ）   d ，使得

A  -1σ，τ +

u，τ - σ，

v + cu，v =

〈f，v〉， （v，τ）∈ H  1  0（ Ω ）×   L  2（ Ω ）   d  （4）

为避免证明LBB条件，我们引入一个非对称的残差项 （ A  -1 σ  h+

u  h，τ-A

v ）  T ，施加在每个单元  T∈ T  h 上  ［12］，并将 （4） 式中的梯度项

替换为弱梯度

w ，得到下面的协调间断有限元数值格式：

格式1  （协调间断有限元格式） 求 （ u  h， σ  h）∈ V  h× W  h ，使得

B  u  h， σ  h;v，τ =〈f，v〉，  v，τ ∈ V  h× W  h  （5）

其中

B  u  h， σ  h;v，τ =   A  -1 σ  h，τ    T  h+

w u  h，τ    T  h-

σ   h，

wv     T   h+  c u   h，v     T   h-∑  T∈ T   h δ（ A  -1 σ   h+

w u   h，τ-A

wv ）   T.

下面证明上述数值格式有存在唯一解. 首先定义如下带有函数 κ=κ（x）>0 形式的内积：

（φ，ψ）   T，κ=∫   Tκ（x）φψ d x，

相应的范数为

φ   T，κ= （φ，φ）  1/2  T，κ .

为证明格式1的解的存在唯一性，先定义与双线性泛函 B  u  h， σ  h;v，τ  对应的稳定化范数    ·，·    δ ：

v，τ     2  δ=∑  T∈ T   h    v    2  T，c+δ

wv    2  T，A+  τ    2   A  -1    （6）

另外，再定义  V  h 中的两个半范数   v   h，1 和   v   h，2 ：

v    2  h，1=∑  T∈ T   h δ

wv    2  T，A，

v    2  h，2=∑  T∈ T   h δ

wv    2  T，A+∑  e∈ F   h δ h  -1   e   v     2  e，A.

v   h，1 为协调间断有限元法下的半范数  ［5］，   v   h，2 为修正的弱有限元方法下的半范数  ［4］. 这两种范数是等价的  ［6］，即存在常数  C  1， C  2 ，使得

C  1  v   h，2≤  v   h，1≤ C  2  v   h，2， v∈ V  h  （7）

下面的引理表明式（5）中的双线性泛函  B（ u  h， σ  h;v，τ） 的强制性.

引理2.1   设 0<δ<1 . 存在常数  c  δ ，使得对任意 （v，τ）∈ V  h× W  h ，下列不等式成立：

B（v，τ;v，τ）≥ c  δ   v，τ    2  δ.

证明  由式（6）中范数    ·，·    δ 的定义直接可得

B v，τ;v，τ =∑  T∈ T   h  1-δ （ A  -1τ，τ ）  T+

∑  T∈ T   h δ（A

wv，

wv ）   T+  cv，v     T   h≥

c   δ   v，τ     2  δ.

证毕.

下面的引理证明了式（5）中双线性泛函  B（ u  h， σ  h;v，τ） 的有界性.

引理2.2   设 0<δ<1.  存在常数  C  δ ，使得对任意   v  1， τ  1 ∈ V  h× W  h ， （ v  2， τ  2）∈ V  h× W  h ，有

B（ v  1， τ  1; v  2， τ  2）≤ C  δ  （ v  1， τ  1）   δ  （ v  2， τ  2）   δ.

证明  由离散的Cauch-Schwarz不等式，有

B  v  1， τ  1; v  2， τ  2 ≤

［4   A  -1 τ  1， τ  1     T   h+

A

w v  1，

w v  1     T   h+2∑  T∈ T   h δ（A

w v  1，

w v  1 ）   T+

c v  1， v  1     T   h ］   1 2 *

［4   A  -1 τ  2， τ  2     T   h+

A

w v  2，

w v  2     T   h+2∑  T∈ T   h δ（A

w v  2，

w v  2 ）   T+

c v  2， v  2     T   h ］   1 2 ≤［4   A  -1 τ  1， τ  1     T   h+

A

w v  1，

w v  1     T   h+C∑  T∈ T   h δ（A

w v  1，

w v  1 ）   T+

c v  1， v  1    T   h ］   1 2 *［4   A  -1 τ  2， τ  2     T   h+

A

w v  2，

w v  2     T   h+ C∑  T∈ T   h δ（A

w v  2，

w v  2 ）  T+  c v  2， v  2     T   h ］   1 2 ≤

C   δ  （ v  1， τ  1）    δ  （ v  2， τ  2）    δ.

证毕.

由引理2.1和2.2可知，双线性泛函  B（ u  h， σ  h;v，τ） 是强制有界的，且 （5） 式右端的泛函 f 是线性连续的. 进而，根据引理2.1、2.2及Lax-Milgram定理  ［13，14］，格式（5）存在唯一解.

2.2 误差分析

为给出式（5）的误差分析，先定义三个  L  2 -局部投影算子用于误差方程推导. 对任意单元 T∈ T  h ，令   Π   j  h ，   Π   k  h 和   Π   k-1  h 分别为  L  2（T） 到 ［ P  j（T） ］  d ，   P  k（T） 及 ［ P  k-1（T） ］  d 上的投影算子，其中算子   Π   j  h 满足如下引理中的方程.

引理2.3   ［11］  若 v∈ H  1  0（ Ω ），  则在单元 T∈ T  h 上下列方程成立.

wv=  Π   j  h

v  （8）

下面我们给出协调间断有限元格式（5）的误差方程. 令

ε  h=  Π   k  hu- u  h， ε  h=  Π   k-1  hσ- σ  h .

则下面引理中的误差方程成立.

引理2.4   设 （u，σ） 是问题（2）的精确解， （ u  h， σ  h）∈ V  h× W  h 是协调间断元格式（5）的解. 则对任意 （v，τ）∈ V  h× W  h ，下面的误差方程成立.

B  ε  h， ε  h，v，τ = l  u，σ v，τ ，

其中

l   u，σ（v，τ）= l  1（u，σ;v，τ）- l  2（u，τ）- l  3（σ，v），

l  1（u，σ;v，τ）=∑  T∈ T   h δ（ A  -1 ε   h-（

u-  Π    j  h

u），

τ-A

wv ）   T，

l  2（u，τ）= 〈u-  Π    k  hu，τ·n〉     T   h，

l  3 σ，v = 〈 σ-  Π    k-1  hσ ·n，v- v 〉     T   h.

证明  考虑弱形式 B u，σ;v，τ = f，v . 下面我们用分布积分公式和弱梯度

wv 的定义推出式（1）和（5）中各对应项的误差方程.

u，τ =- u，

·τ + 〈u，τ·n〉    T  h=

-   Π   k  hu，

·τ + 〈u，τ·n〉    T  h=

（

w  Π   k  hu，τ）+ l  2（u，τ），

·σ，v =- σ，

v + 〈σ·n，v〉    T  h=

-   Π   k-1  hσ，

v + 〈σ·n，v- v 〉    T  h=

Π   k-1  hσ，v - 〈  Π   k-1  hσ·n，v〉    T  h+

〈σ·n，v- v 〉    T  h=-   Π   k-1  hσ，

w·v -

〈  Π   k-1  hσ·n，v- v 〉    T  h+ 〈σ·n，v- v 〉    T  h=

-（  Π   k-1  hσ，

w·v）+ l  3（σ，v），

∑  T∈ T   h δ（ A  -1  Π    k-1  hσ+

w  Π    k  hu，τ-A

wv ）   T=

∑  T∈ T   h δ（ A  -1  Π    k-1  hσ-σ-（

u-  Π    j  h

u），τ-

A

wv ）   T= l  1（u，σ;v，τ）.

将   Π   k  hu 和   Π   k-1  hσ 代入式（1），  u  h 和  σ  h 代入式（5），联立可得

B  ε  h， ε  h，v，τ = l  1 u，σ;v，τ - l  2 u，τ -

l  3（σ，v）.

下面我们给出精确解 （u，σ） 和离散解 （ u  h， σ  h） 在范数   （·，·）   δ 下的最优收敛阶误差估计. 本节先给出迹不等式和跳量  v  和均值  v  的关系，用于下面的误差分析.

引理2.5    ［2］ 对任意函数 φ∈ H  1（T）， 下面的迹不等式成立.

φ   2  e≤C  h  -1  T  φ   2  T+ h  T

φ   2  T   （9）

对于任意 v∈ V  h ， v 的跳量  v  和均值  v  满足如下关系.

v- v    e=  v   e， e?  Ω ，

v-{v}   e= 1 2    v    e， e∈ F  0.

下面我们给出引理2.5中式  l  1 u，σ;v，τ ，  l  2（u，τ） 和  l  3 σ，v  的误差估计.

引理2.6   若反应系数 c 在每个 单元 T 上是一次可微的连续函数，那么对于 （u，σ）∈ H  k+1（ Ω ）× ［ H  k（ Ω ） ］  d ，下列估计成立.

| l  1（u，σ;v，τ）|≤C h  k（  σ   k+    u   k+1）  （v，τ）   δ，

| l  2（u，τ）|≤C h  k  u   k+1  （v，τ）   δ，

l  3 σ，v  ≤C h  k  σ   k   v，τ    δ.

证明  由Cauchy-Schwarz不等式和迹不等式（9），有

l  1 u，σ;v，τ  =|∑  T∈ T   h δ（ A  -1 ε   h-（

u-

Π    j  h

u），τ-A

wv ）   T|≤C（∑  T∈ T   h δ（   ε   h   T+

u-  Π     j  h

u    T））   1 2 ×  ∑  T∈ T   h   τ-A

wv    T    1 2 ≤

C h   k（  σ    k+  u    k+1）  （v，τ）    δ.

同理，

l  2 u，τ  = ∑  T∈ T   h  〈u-  Π    k  hu，τ·n〉    T ≤

C  ∑  T∈ T   h  h   -1  T  u-  Π    k  hu    2   T    1 2 ×

∑  T∈ T   h  h   T  τ    2   T   1/2≤C h   k  u    k+1  （v，τ）    δ.

由Cauchy-Schwarz不等式，迹不等式（9）和式（7）有

l  3 σ，v  =

|∑  T∈ T   h  〈 σ-  Π    k-1  hσ ·n，v- v 〉    T|≤

C  ∑  T∈ T   h  h   T  σ-  Π    k-1  hσ   2    T   1/2×

∑  T∈ T   h  h  -1  T   v     2  e   1/2≤

C h  k  σ    k  （v，τ）   δ.

定理2.7   设 （u，σ） 是式（2）的精确解， （ u  h， σ  h）∈   V  h× W  h 是协调间断有限元格式（3）的解，则下列误差估计不等式成立：

Π   k  hu- u  h，  Π   k-1  hσ- σ  h    δ≤

C h  k（  u   k+1+  σ   k）  （10）

证明  令 v= ε  h ，  τ= ε  h . 我们有

B  δ（  h，  h;  h，  h）= l  u，σ（v，τ），

其中

l  u，σ（v，τ）= l  1（u，σ;  h，  h）- l  2（u，  h）- l  3（σ，  h）.

由引理2.1和2.2，下列不等式成立：

（  h，  h）   δ≤C h  k（  u   k+1+  σ   k），

从而推出式（10）. 证毕.

定理2.8   设 （u，σ） 是式（2）的精确解， （ u  h， σ  h）∈   V  h× W  h 是协调间断有限元格式（3）的解，则下列误差估计不等式成立：

u- u  h，σ- σ  h    δ≤C h  k   u   k+1+  σ   k   （11）

证明  由定理2.3和  L  2 投影算子   Π   k  h 和   Π   k-1  h 的逼近性质可直接得.

2.3 数值算例

对二维定常扩散反应方程，我们用两个算例验证格式1的有效性. 针对以下两个算例，我们对区域Ω=［0，1］×［0，1］采用一致三角形剖分.

例2.9   在式（1）中令 A = I，I 为单位矩阵， c=1 ， σ=-A

u ，则有

u= sin （π x  1） sin （π x  2），

σ=-π  cos （π x  1） sin （π x  2）  sin （π x  1） cos （π x  2） .

取 k=1可 得到下表1中的误差.

取 h=1/16，1/32，1/64，1/128 ，我们得到 A= I ，c=1 时定常扩散反应方程在能量范数   （·，·）   δ 下的误差和收敛阶. 表1显示，原始变量 u 采用P  1元，通量 σ 采用P  0元，数值结果的收敛阶为一阶，验证了数值格式的有效性.

例2.10   在式（1）中令 A = 1e-3  I，I 为单位矩阵， c=1 ，则方程变为反应占优的扩散反应方程.

取 k=1 可得到下表2中的误差.

表2显示，原始变量 u 采用P  1元，通量 σ 采用P  0元， h=1/64，1/128 时，收敛阶可以达到一阶，证明数值格式的鲁棒性.

3 非定常扩散反应方程情形

3.1 全离散协调间断有限元格式

本节中我们将给出非定常扩散反应扩散方程的形式和全离散协调间断有限元格式.

考虑如下问题：求函数 u=u（x，t） ，满足方程

u  t-

· A

u +cu=f，（x，t）∈ Ω ×（0， Γ ］，

u=0， x，t ∈  Ω × 0， Γ  ，

u x，0 = u  0 x ，x∈ Ω       （12）

其中  Γ <+∞ 为最终时刻，  u  0（x） 为 u（x，t） 的初始值.

引入通量 σ=-A

u ，则方程（12）可以化成以下等价的混合形式：

u+ A  -1σ=0，（x，t）∈ Ω ×（0， Γ ］，

u  t+

·σ+cu=f， x，t ∈  Ω × 0， Γ  ，

u x，0 = u  0 x ，x∈ Ω       （13）

令 Δt 为时间步长，  t  i=iΔt，i=1，…，N ，  u  i  h= u  h（ t  i） ，  σ  i  h= σ  h（ t  i） 和  f  i=f（·， t  i） . 在 t= t  i 时，本文采用向后欧拉差分格式，则可得到下面的全离散协调间断有限元数值格式：

格式2  （全离散协调间断有限元格式） 求 （ u  i  h， σ  i  h）∈ V  h× W  h，i=1，2，…，N ，对任意  v，τ ∈ V  h× W  h ，使得

t  u  i  h，v    T  h+B  u  i  h， σ  i  h;v，τ =〈 f  i，v〉，

u  0  h= Q  h u  0  （14）

或

u  i  h，v    T  h+ΔtB  u  i  h， σ  i  h;v，τ =

Δt〈 f  i，v〉+   u  i-1  h，v    T  h， u  0  h= Q  h u  0，

其中

B  u   i  h， σ   i  h;v，τ =   A  -1 σ   i  h，τ     T   h+

w u   h，τ     T   h-

σ   i  h，

wv     T   h+  c u   i  h，v     T   h-∑  T∈ T   h δ（ A  -1 σ   i  h+

w u   i  h，τ-A

wv ）   T.

3.2 误差分析

为便于分析误差，先定义两个椭圆投影.   R  k  h 为  L  2（0， Γ ; H  1  0（ Ω ）） 到  V  h 上的椭圆投影和  R  k-1  h 为  L  2（0， Γ ;［ L  2（ Ω ） ］  d） 到  W  h 上的椭圆投影，这两个椭圆投影类似于文献［15］中的Wheeler投影.

设 （u，σ）∈ L  2（0， Γ ; H  1  0（ Ω ））× L  2（0， Γ ;［ L  2（ Ω ） ］  d） 为问题（13）的精确解. 对任意 t∈（0， Γ ］ ，我们考虑如下问题：

求   R  k  hu， R  k-1  hσ ∈ V  h× W  h ，对任意 （v，τ）∈    V  h× W  h ，使得

B  R  k  hu， R  k-1  hσ;v，τ =B u，σ;v，τ   （15）

式（15）可以视为定常扩散反应方程协调间断元方法的数值格式（5），所以 （ R  k  hu， R  k-1  hσ） 由式（15）唯一解出. 因此，由引理2.1可以立即推出如下的引理.

引理3.1   设 （u，σ）∈ L  2（0， Γ ; H  1  0（ Ω ）∩   H  k+1（ Ω ））× L  2（0， Γ ;［ H  k（ Ω ） ］  d） 为问题（13）的精确解， （ R  k  hu， R  k-1  hσ）∈ V  h× W  h 为问题（15）的解，则存在常数 C ，使得下列不等式成立.

（  Π   k  hu- R  k  hu，  Π   k-1  hσ- R  k-1  hσ）   δ≤

C h  k（  u   k+1+  σ   k）.

若进一步假设  u  t∈ H  1  0（ Ω ）∩ H  k+1（ Ω ） ，则有

Π   k  h u  t- R  k  h u  t ≤C h  k+1（   u  t   k+1+  σ   k）.

结合引理3.1和引理2.2，并由三角不等式，我们有如下引理成立.

引理3.2   设 （u，σ） 为问题（13）的精确解，  （ R  k  hu， R  k-1  hσ）∈ V  h× W  h 为问题（15）的解，则对任意 t∈（0， Γ  ， （u（x），σ（x））∈ H  1  0（ Ω ）? H  k+1（ Ω ）× ［ H  k（ Ω ） ］  d ，有如下估计成立.

（u- R  k  hu，σ- R  k-1  hσ）   δ≤

C h  k（  u   k+1+  σ   k）.

下面我们推导全离散协调间断有限元格式的误差估计.令  μ  i= R  k  h u  i- u  i  h ，  ν  i= R  k-1  h σ  i- σ  i  h ， i=1，2，…，N ，则有以下引理成立.

引理3.3   设 （ u  i  h， σ  i  h）∈ V  h× W  h，i=1，2，…，N 为全离散协调间断有限元格式（14）的解， （u，σ）∈ L  2（0， Γ ; H  1  0（ Ω ））× L  2（0， Γ ;［ L  2（ Ω ） ］  d） 为问题（13）的精确解，则对任意  （v，τ）∈ V  h× W  h ，下列等式成立.

t  μ  i，v    T  h+B  μ  i， ν  i;v，τ =-       t  η  i，v    T  h-

（     t  ρ  i，v）   T  h- （ u  i  t-     t  u  i，v）   T  h  （16）

其中  η  i=  Π   k  h u  i- R  k  h u  i 和  ρ  i= u  i-  Π   k  h u  i  h.

证明  由式（12），（14）和（15），有

t  μ  i，v    T  h+B  μ  i， ν  i;v，τ =

t   R  k  h u  i- u  i  h ，v    T  h+

B  R  k  h u  i- u  i  h， R  k-1  h σ  i- σ  i  h;v，τ =

t  R  k  h u  i  h，v    T  h+B  R  k  h u  i， R  k-1  h σ  i;v，τ -〈 f  i，v〉=

t  R  k  h u  i  h，v    T  h+B  u  i， σ  i;v，τ -〈 f  i，v〉=

t  R  k  h u  i  h，v    T  h-   u  i  t，v    T  h=

t  R  k  h u  i  h-     t   Π   k  h u  i  h+     t   Π   k  h u  i  h-     t  u  i+     t  u  i- u  i  t，v    T  h=

- （     t  η  i，v）   T  h- （     t  ρ  i，v）   T  h- （ u  i  t-     t  u  i，v）   T  h.

证毕.

下面的引理描述了  u  i  h- R  k  h u  i 和  σ  i  h- R  k-1  h σ  i，i=1，2，…，N 的最优收敛阶误差估计.

引理3.4   设 （ u  i  h， σ  i  h）∈ V  h× W  h，i=1，2，…，N  为全离散协调间断有限元格式（14）的解， （u，σ）∈  L  2（0， Γ ; H  1  0（ Ω ）∩ H  k+1（ Ω ））× L  2（0， Γ ;   H  k（ Ω ）   d） 为问题（13）的精确解，则对于任意 （v，τ）∈ V  h× W  h ，下列不等式成立.

（ u   i  h- R   k  h u   i， σ   i  h- R   k-1  h σ   i）    δ≤

（ u   0  h- R   k  h u  0， σ   i  h- R   k-1  h σ   i）    δ+

C  h   k+1 ∫    t   i  0（   u   t    k+1+  σ    k） d t+Δt ∫    t   i  0  u   tt  d t    （17）

t （ u   i  h- R   k  h u   i） ≤

u   0  h- R   k  h u  0， σ   i  h- R   k-1  h σ   i    δ+

C  h   k+1 ∫    t   i  0（   u   t    k+1+  σ    k） d t+Δt ∫    t   i  0  u   tt  d t    （18）

证明  在式（16）中令 v=     t  μ  i 和 τ= ν  i 得

t  μ  i   2+B      t  μ  i， ν  i;     t  μ  i， ν  i =

-       t  η  i，     t  μ  i    T  h- （     t  ρ  i，     t  μ  i）   T  h-

（ u  i  t-     t  u  i，     t  μ  i）   T  h   （19）

由式（5）中 B 的强制性、有界性、向后欧拉差分格式和Young不等式，有

B  μ  i， ν  i;     t  μ  i， ν  i =B  μ  i， ν  i;  μ  i- μ  i-1 Δt ， ν  i ≥

1 2Δt   （ μ  i， ν  i）   2  δ- 1 2Δt   （ μ  i-1， ν  i）   2  δ  （20）

对式（19）应用Cauchy-Schwarz 不等式，有

（     t  η  i，     t  μ  i）   T  h≤       t  η  i   2+1/4       t  μ  i   2  （21）

其中

t  η   i = 1 Δt   ∫    t   i   t   i-1  Π    k  h u   t- R   k  h u   t d t ≤

1 Δt  C  1 h   k+1∫   t   i   t   i-1 （   u   t    k+1+  σ    k） d t  （22）

（     t  ρ  i，     t  μ  i）   T  h≤       t  ρ  i   2+1/4       t  μ  i   2  （23）

其中

t  ρ   i = 1 Δt   ∫    t   i   t   i-1  Π    k  h u   t- u   t d t ≤

1 Δt  C  1 h   k+1∫   t   i   t   i-1 （   u   t    k+1+  σ    k） d t  （24）

（ u  i  t-     t  u  i，     t  μ  i）   T  h≤ 1 2    u  i  t-     t  u  i   2+ 1 2        t  μ  i   2  （25）

其中

u   i  t-      t  u   i = 1 Δt   ∫    t   i   t   i-1 t- t   i-1  u   tt（t） d t ≤

∫   t   i   t   i-1   u   tt  d t  （26）

结合上述不等式， 有

μ  i， ν  i    2  δ≤    μ  i-1， ν  i    2  δ+2Δt       t  η  i   2+

2Δt       t  ρ  i   2+Δt   u  i  t-     t  u  i   2.

由此可得

μ  i， ν  i    2  δ≤    μ  i-1， ν  i    2  δ+2Δt       t  η  i   2+

2Δt       t  ρ  i   2+Δt   u  i  t-     t  u  i   2≤

μ   i-2， ν   i-1    2+2Δt ∑  i  j=i-1         t  η   j   2+

∑  i  j=i-1         t  ρ   j   2 +Δt∑  i  j=i-1    u   j  t-      t  u   j   2≤

≤   μ  0   2+2Δt∑  i  j=1         t  η   j   2+

2Δt∑  i  j=1         t  ρ   j   2+Δt∑  i  j=1    u   j  t-      t  u   j   2.

结合式（19）～（26）可以推出式（17）.

同理，由Young不等式有

（ u  i  t-     t  u  i，     t  μ  i）   T  h≤ ‖ u  i  t-     t  u  i‖  2+

1/4 ‖     t  μ  i‖  2  （27）

联立式（19）～（27）可得式（18）. 证毕.

由引理3.1和引理3.3及三角不等式可以推出下面的定理.

定理3.4   设 （ u  i  h， σ  i  h）∈ V  h× W  h，i=1，2，…，N  为全离散协 调间断有限元格式（14）的解， （u，σ）∈  L  2（0， Γ ; H  1  0（ Ω ）∩ H  k+1（ Ω ））× L  2（0， Γ ;   H  k（ Ω ）   d） 为式（13）的精确解，则下列估计成立.

t （ u   i- u   i  h） +  （ u   i- u   i  h， σ   i- σ   i  h）    δ≤

（ u   0  h- R   k  h u  0，0）    δ+

C  h   k（  u    k+1+  σ    k）+Δt ∫    t   i  0  u   tt  d t .

3.3 数值算例

例3.5   在式（11）中令 A = I，I 为单位矩阵，  c=1 ，则有

u=16t x  1（1- x  1） x  2（1- x  2），

σ=16t  （2 x  1-1） x  2（1- x  2）

（2 x  2-1） x  1（1- x  1）  .

取 k=1 可以得到下表3中的误差.

分别取 Δt=1/16，1/32，1/64 ， h=1/16，1/32，1/64 ，可以得到 当A= I ，c=1 时非定常定常扩散反应方程在能量范数   （·，·）   δ 下的误差和收敛阶. 从表3可以看出，原始变量 u 采用P  1元，通量 σ 采用P  0元，数值结果的收敛阶为一阶，数值格式有效.

例3.6   令 A =（1e-3） I，I 为单位矩阵， c=1 ，则式（12）变为反应占优的扩散反应方程.

取 k=1 可以得到下表4中的误差.

从表4可以看出，原始变量 u 采用P  1元，通量 σ 采用P  0元， Δt=1/32，1/64，h=1/32，1/64 时，收敛阶可以达到一阶，因而本节中的数值格式具有鲁棒性.

4 结 论

本文对定常和非定常的反应扩散方程提出了一个全离散协调间断有限元格式. 算例表明该格式有效且鲁棒.

值得注意的是，协调间断有限元法虽然形式简洁，但在计算弱梯度时因使用了更高次数的弱梯度有限元空间而可能导致计算弱梯度的系数矩阵条件数很大，病态较为严重，从而使得求解较为困难. 此外，可令 δ=0 且去掉格式中的非对称稳定项，进而通过验证LBB条件来证明去掉非对称稳定项后的格式也有唯一解，且形式更简单.

参考文献：

［1］   Wang J， Ye X. A weak Galerkin finite element method for second-order elliptic problems ［J］. J Comput Appl Math， 2013， 241： 103.

［2］  Wang J， Ye X. A weak Galerkin mixed finite element method for second order elliptic problems ［J］. Math Comput， 2014， 83： 2101.

［3］  Mu L， Wang J， Ye X. A least-squares-based weak Galerkin finite element method for second order elliptic equations ［J］. SIAM J Sci Comput， 2017， 39： A1531.

［4］  Wang X， Malluwawadu N S， Gao F，  et al . A modified weak Galerkin finite element method ［J］. J Comput Appl Math， 2014， 271： 319.

［5］  Ye X， Zhang S. A conforming discontinuous Galerkin finite element method ［J］. Int J Numer Anal Model， 2020， 17： 110.

［6］  Ye X， Zhang S. A conforming discontinuous Galerkin finite element method （II） ［J］. Int J Numer Anal Model， 2020， 17： 281.

［7］  Ye X， Zhang S. A conforming discontinuous Galerkin finite element method （III） ［J］. Int J Numer Anal Model， 2020， 17： 794.

［8］  Mu L， Ye X， Zhang S. Development of pressure-robust discontinuous Galerkin finite element methods for the Stokes problem ［J］. J Sci Comput， 2021， 89： 1.

［9］  Pehlivanov A I， Carey G F， Lazarov R D. Least-squares mixed finite elements for second-order elliptic problems ［J］. SIAM J Numer Anal， 1994， 31： 1368.

［10］  李西， 罗家福， 冯民富. 非定常Navier-Stokes方程的一种非线性局部投影稳定化有限元方法［J］. 四川大学学报：自然科学版， 2021， 58： 031002.

［11］ Ye X， Zhang S. A conforming discontinuous Galerkin finite element method for the Stokes problem onpolytopal meshes ［J］. Int J Numer Meth Fluids， 2021， 93： 1913.

［12］ Fu H， Guo H， Hou J，  et al . A stabilized mixed finite element method for steady and unsteady reaction-diffusion equations ［J］. Comput Meth Appl Mech E， 2016， 304： 102.

［13］ 李开泰，黄艾香，黄庆怀. 有限元方法及其应用［M］. 西安： 西安交通大学出版社， 1988.

［14］ 张世清. 泛函分析及其应用［M］. 北京： 科学出版社， 2012.

［15］ Wang X，Zhai Q， Zhang R. The weak Galerkin method for solving the incompressible Brinkman flow ［J］. J Comput Appl Math， 2016， 307： 13.