顧子兵 胡朝浪 楊榮奎 馮民富

摘 要 ???：對于定常擴散反應(yīng)方程的原始混合變分形式，本文基于間斷有限元法和最小二乘法思想給出了一種新的協(xié)調(diào)間斷有限元格式，并將其推廣應(yīng)用于非定常擴散反應(yīng)方程，進(jìn)而建立了全離散協(xié)調(diào)間斷有限元格式. 該格式不僅能避免LBB條件，而且適用于多邊形網(wǎng)格. 在能量范數(shù)下，本文得到了原始變量和通量的最優(yōu)誤差估計.最后，針對定常和非定常的擴散反應(yīng)方程，本文分別以一般情形和對流占優(yōu)情形下的數(shù)值算例驗證了方法的有效性.

關(guān)鍵詞 ：協(xié)調(diào)間斷有限元； 原始混合形式； 最小二乘法

中圖分類號 ： O241.82 文獻(xiàn)標(biāo)識碼 ：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 引 言

本文考慮如下擴散反應(yīng)問題：求函數(shù) u=u（x） ，使得

-??（1）

其中，Ω是 ?Euclid Math TwoRA@

d d=2，3 ?中的一個多邊形區(qū)域，其Lipschitz邊界為 ??Ω， A 為正定矩陣， c≥0 . 該問題常被用于描述溫度、壓力及化學(xué)物質(zhì)的濃度等物理量的演化.

有限元法是求解問題（1）的一種有效手段. 用經(jīng)典的有限元法可以直接求得問題（1）的離散解. 在實際問題中， u 的梯度

u 往往有重要的物理意義，如在流動問題中

u 表示流體的速度. 為求 u 和

u， 2013年Wang和Ye ?［1］對二階橢圓問題提出了弱有限元法. 該方法區(qū)別于間斷有限元法的地方在于它將一個不連續(xù)的函數(shù) v 在單元上的取值分別為內(nèi)部 ?v ?0 和邊界 ?v ?b ，在構(gòu)造有限元時分別在內(nèi)部和邊界上構(gòu)造函數(shù)，用弱梯度

w· 代替梯度

. 弱有限元法具有和間斷有限元法相似的優(yōu)點，且形函數(shù)選取及網(wǎng)格剖分多樣. 2014年，Wang和Ye ?［2］針對帶Dirichlet邊界的擴散方程的對偶混合形式提出了弱Galerkin混合有限元方法.該方法對主要變量和通量都有非常精確地近似，且可以應(yīng)用于多邊形和多面體網(wǎng)格. 2017年，Mu等 ?［3］提出了最小二乘的弱有限元方法，該方法可以被用于奇異攝動擴散反應(yīng)方程和各向異性、不均勻性多孔介質(zhì)中流體的流動問題，具有較好的魯棒性和靈活性. 值得注意的是，上述工作都是基于方程的對偶混合變分形式提出的方法.

2014年，Wang等 ?［4］提出了修正的弱有限元法.該方法將弱有限元方法中的邊界函數(shù) ?v ?b 替換為 ?v ?，通過單元內(nèi)部函數(shù) v 確定邊界函數(shù) ?v ?， 減少了邊界自由度，提高了計算效率. 2020年，Ye等 ?［5-8］提出了協(xié)調(diào)間斷元方法.該方法屬于弱有限元方法，邊界函數(shù)仍然為 ?v ?. 此外，該方法還具備協(xié)調(diào)有限元和間斷有限元方法的特征，形式簡潔，單元形函數(shù)采用間斷函數(shù)，且可以通過調(diào)整有限元空間將該方法應(yīng)用于多邊形和多面體網(wǎng)格.

基于上述工作，本文考慮方程（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），本文采用協(xié)調(diào)間斷有限元方法結(jié)合最小二乘法 ?［3，9］，給出了一種新的協(xié)調(diào)間斷有限元格式. 通過引入最小二乘的殘差項，該方法不但避免了LBB條件 ?［10］的限制， 而且可被應(yīng)用于多邊形網(wǎng)格.

本文結(jié)構(gòu)如下. 第2節(jié)針對定常擴散反應(yīng)方程給出了一種協(xié)調(diào)間斷有限元方法的數(shù)值格式，并給出理論分析. 在此基礎(chǔ)上，本文也給出一般情況和反應(yīng)占優(yōu)時的數(shù)值算例，驗證了該格式的有效性. 在第3節(jié)中，本文將第2節(jié)中的協(xié)調(diào)間斷有限元格式推廣應(yīng)用于非定常擴散反應(yīng)方程，給出了一種全離散協(xié)調(diào)間斷有限元格式，并通過一般情況和反應(yīng)占優(yōu)時的數(shù)值算例驗證了該格式的有效性. 第4節(jié)總結(jié)本文所得結(jié)果.

2 定常擴散反應(yīng)方程情形

2.1 協(xié)調(diào)間斷有限元格式

我們以 ?W ?m，p（ Ω ） ， ?W ?m，p ?0（ Ω ） 表示Ω上的 m -階Sobolev空間. 特別地，當(dāng) p =2時，

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

·，· ??m ， ??· ??m ， ??· ??m 分別為 ?H ?m（ Ω ） 的內(nèi)積、范數(shù)和半范數(shù).

令 ?T ?h 為區(qū)域 ?Ω ?的有限元剖分，該剖分由凸且正則的多邊形或多面體組成. 記 ?F ?h 為 ?T ?h 的所有邊或面的集合， ?F ?0 ?h= F ?h＼ ?Ω ?表示所有內(nèi)邊的集合. 為了簡便，分別定義 ?T ?h 和 ??T ?h 上的內(nèi)積為

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.

為了簡便表示函數(shù)的弱梯度，用 v 而不用 ?v ?h 表示有限元空間 ?V ?h 的任意函數(shù).

對每個單元 T∈ T ?h ，記 ?h ?T 為其直徑， ?T ?h 的網(wǎng)格大小為 h= ?max ??T∈ T ?h ???h ?T . 定義弱函數(shù) v ，

v= v， 在單元T內(nèi)部，

{v}， 在邊界 T上.

設(shè)整數(shù) k≥1 ，記 ?P ?k（T） 為定義在 T 上次數(shù)不超過 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.

對于函數(shù) 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 的選取依賴于多邊形的邊數(shù). 對于三角形網(wǎng)格， j 可以取 k+1 .通常情況下， j=n+k-1 ，其中 n 為多邊形的邊數(shù) ?［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 ，得到下面的協(xié)調(diào)間斷有限元數(shù)值格式：

格式1 ?（協(xié)調(diào)間斷有限元格式） 求 （ 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.

下面證明上述數(shù)值格式有存在唯一解. 首先定義如下帶有函數(shù) κ=κ（x）>0 形式的內(nèi)積：

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

相應(yīng)的范數(shù)為

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

為證明格式1的解的存在唯一性，先定義與雙線性泛函 B ?u ?h， σ ?h;v，τ ?對應(yīng)的穩(wěn)定化范數(shù) ???·，· ???δ ：

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

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

另外，再定義 ?V ?h 中的兩個半范數(shù) ??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 為協(xié)調(diào)間斷有限元法下的半范數(shù) ?［5］， ??v ??h，2 為修正的弱有限元方法下的半范數(shù) ?［4］. 這兩種范數(shù)是等價的 ?［6］，即存在常數(shù) ?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 ??設(shè) 0<δ<1 . 存在常數(shù) ?c ?δ ，使得對任意 （v，τ）∈ V ?h× W ?h ，下列不等式成立：

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

證明 ?由式（6）中范數(shù) ???·，· ???δ 的定義直接可得

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 ??設(shè) 0<δ<1. ?存在常數(shù) ?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 是線性連續(xù)的. 進(jìn)而，根據(jù)引理2.1、2.2及Lax-Milgram定理 ?［13，14］，格式（5）存在唯一解.

2.2 誤差分析

為給出式（5）的誤差分析，先定義三個 ?L ?2 -局部投影算子用于誤差方程推導(dǎo). 對任意單元 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）

下面我們給出協(xié)調(diào)間斷有限元格式（5）的誤差方程. 令

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

則下面引理中的誤差方程成立.

引理2.4 ??設(shè) （u，σ） 是問題（2）的精確解， （ u ?h， σ ?h）∈ V ?h× W ?h 是協(xié)調(diào)間斷元格式（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）中各對應(yīng)項的誤差方程.

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），聯(lián)立可得

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

l ?3（σ，v）.

下面我們給出精確解 （u，σ） 和離散解 （ u ?h， σ ?h） 在范數(shù) ??（·，·） ??δ 下的最優(yōu)收斂階誤差估計. 本節(jié)先給出跡不等式和跳量 ?v ?和均值 ?v ?的關(guān)系，用于下面的誤差分析.

引理2.5 ???［2］ 對任意函數(shù) φ∈ H ?1（T）， 下面的跡不等式成立.

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

φ ??2 ?T ??（9）

對于任意 v∈ V ?h ， v 的跳量 ?v ?和均值 ?v ?滿足如下關(guān)系.

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 ??若反應(yīng)系數(shù) c 在每個 單元 T 上是一次可微的連續(xù)函數(shù)，那么對于 （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 ??設(shè) （u，σ） 是式（2）的精確解， （ u ?h， σ ?h）∈ ??V ?h× W ?h 是協(xié)調(diào)間斷有限元格式（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 ??設(shè) （u，σ） 是式（2）的精確解， （ u ?h， σ ?h）∈ ??V ?h× W ?h 是協(xié)調(diào)間斷有限元格式（3）的解，則下列誤差估計不等式成立：

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

證明 ?由定理2.3和 ?L ?2 投影算子 ??Π ??k ?h 和 ??Π ??k-1 ?h 的逼近性質(zhì)可直接得.

2.3 數(shù)值算例

對二維定常擴散反應(yīng)方程，我們用兩個算例驗證格式1的有效性. 針對以下兩個算例，我們對區(qū)域Ω=［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 時定常擴散反應(yīng)方程在能量范數(shù) ??（·，·） ??δ 下的誤差和收斂階. 表1顯示，原始變量 u 采用P ?1元，通量 σ 采用P ?0元，數(shù)值結(jié)果的收斂階為一階，驗證了數(shù)值格式的有效性.

例2.10 ??在式（1）中令 A = 1e-3 ?I，I 為單位矩陣， c=1 ，則方程變?yōu)榉磻?yīng)占優(yōu)的擴散反應(yīng)方程.

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

表2顯示，原始變量 u 采用P ?1元，通量 σ 采用P ?0元， h=1/64，1/128 時，收斂階可以達(dá)到一階，證明數(shù)值格式的魯棒性.

3 非定常擴散反應(yīng)方程情形

3.1 全離散協(xié)調(diào)間斷有限元格式

本節(jié)中我們將給出非定常擴散反應(yīng)擴散方程的形式和全離散協(xié)調(diào)間斷有限元格式.

考慮如下問題：求函數(shù) 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 時，本文采用向后歐拉差分格式，則可得到下面的全離散協(xié)調(diào)間斷有限元數(shù)值格式：

格式2 ?（全離散協(xié)調(diào)間斷有限元格式） 求 （ 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 上的橢圓投影，這兩個橢圓投影類似于文獻(xiàn)［15］中的Wheeler投影.

設(shè) （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）可以視為定常擴散反應(yīng)方程協(xié)調(diào)間斷元方法的數(shù)值格式（5），所以 （ R ?k ?hu， R ?k-1 ?hσ） 由式（15）唯一解出. 因此，由引理2.1可以立即推出如下的引理.

引理3.1 ??設(shè) （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）的解，則存在常數(shù) C ，使得下列不等式成立.

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

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

若進(jìn)一步假設(shè) ?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）.

結(jié)合引理3.1和引理2.2，并由三角不等式，我們有如下引理成立.

引理3.2 ??設(shè) （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）.

下面我們推導(dǎo)全離散協(xié)調(diào)間斷有限元格式的誤差估計.令 ?μ ?i= R ?k ?h u ?i- u ?i ?h ， ?ν ?i= R ?k-1 ?h σ ?i- σ ?i ?h ， i=1，2，…，N ，則有以下引理成立.

引理3.3 ??設(shè) （ u ?i ?h， σ ?i ?h）∈ V ?h× W ?h，i=1，2，…，N 為全離散協(xié)調(diào)間斷有限元格式（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 的最優(yōu)收斂階誤差估計.

引理3.4 ??設(shè) （ u ?i ?h， σ ?i ?h）∈ V ?h× W ?h，i=1，2，…，N ?為全離散協(xié)調(diào)間斷有限元格式（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）應(yīng)用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）

結(jié)合上述不等式， 有

μ ?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.

結(jié)合式（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）

聯(lián)立式（19）～（27）可得式（18）. 證畢.

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

定理3.4 ??設(shè) （ u ?i ?h， σ ?i ?h）∈ V ?h× W ?h，i=1，2，…，N ?為全離散協(xié) 調(diào)間斷有限元格式（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 數(shù)值算例

例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 ，可以得到 當(dāng)A= I ，c=1 時非定常定常擴散反應(yīng)方程在能量范數(shù) ??（·，·） ??δ 下的誤差和收斂階. 從表3可以看出，原始變量 u 采用P ?1元，通量 σ 采用P ?0元，數(shù)值結(jié)果的收斂階為一階，數(shù)值格式有效.

例3.6 ??令 A =（1e-3） I，I 為單位矩陣， c=1 ，則式（12）變?yōu)榉磻?yīng)占優(yōu)的擴散反應(yīng)方程.

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

從表4可以看出，原始變量 u 采用P ?1元，通量 σ 采用P ?0元， Δt=1/32，1/64，h=1/32，1/64 時，收斂階可以達(dá)到一階，因而本節(jié)中的數(shù)值格式具有魯棒性.

4 結(jié) 論

本文對定常和非定常的反應(yīng)擴散方程提出了一個全離散協(xié)調(diào)間斷有限元格式. 算例表明該格式有效且魯棒.

值得注意的是，協(xié)調(diào)間斷有限元法雖然形式簡潔，但在計算弱梯度時因使用了更高次數(shù)的弱梯度有限元空間而可能導(dǎo)致計算弱梯度的系數(shù)矩陣條件數(shù)很大，病態(tài)較為嚴(yán)重，從而使得求解較為困難. 此外，可令 δ=0 且去掉格式中的非對稱穩(wěn)定項，進(jìn)而通過驗證LBB條件來證明去掉非對稱穩(wěn)定項后的格式也有唯一解，且形式更簡單.

參考文獻(xiàn)：

［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方程的一種非線性局部投影穩(wěn)定化有限元方法［J］. 四川大學(xué)學(xué)報：自然科學(xué)版， 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］ 李開泰，黃艾香，黃慶懷. 有限元方法及其應(yīng)用［M］. 西安： 西安交通大學(xué)出版社， 1988.

［14］ 張世清. 泛函分析及其應(yīng)用［M］. 北京： 科學(xué)出版社， 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.