劉軒宇,羅 鯤,王 皓
(四川大學(xué)數(shù)學(xué)學(xué)院,成都 610064)
考慮下面的拋物型積分微分方程:
其中Ω?Rd(d=2,3) 是一個(gè)有界凸多邊形(多面體),其邊界記為 ?Ω.在物理和工程領(lǐng)域,同時(shí)帶時(shí)間變量和空間變量的系統(tǒng)非常常見.在研究其中一些問題時(shí),必須考慮之前的狀態(tài)對現(xiàn)在狀態(tài)的影響.拋物型積分微分方程應(yīng)運(yùn)而生.這類方程在熱傳導(dǎo)、核反應(yīng)堆動力學(xué)、熱彈性力學(xué)等領(lǐng)域有廣泛的應(yīng)用[1-24].針對這類方程,近幾十年來出現(xiàn)了很多數(shù)值方法,如混合有限元方法[7,13,18],譜方法[8],全離散自適應(yīng)有限元方法[17],間斷Galerkin方法[14,15],最小二乘Galerkin方法[9],緊致的交替隱式差分格式[16]等.
弱Galerkin (Weak Galerkin,簡稱WG)有限元法是近幾年發(fā)展起來的用于求解偏微分方程的一類有限元方法[25-29],最早由Wang和Ye針對二階橢圓方程提出[21,23].WG方法通過引入對間斷函數(shù)適用的弱微分(弱梯度,弱散度)算子,使相應(yīng)變分方程容許間斷的逼近函數(shù).此方法既有間斷有限元法網(wǎng)格剖分靈活的特點(diǎn),也有局部消去的特性,減少了變量個(gè)數(shù).近年來,WG方法更是被推廣應(yīng)用到多種偏微分方程的數(shù)值求解,如線彈性問題[3],Stokes方程[2,20,22,25,27],Ossen方程[12],Biot固結(jié)問題[5],對流擴(kuò)散反應(yīng)方程[1],重調(diào)和方程[1]等.
針對二維拋物型積分微分方程,文獻(xiàn)[29]提出了一類WG有限元方法,給出了半離散格式與全離散格式的解的存在唯一性,導(dǎo)出了相應(yīng)的誤差估計(jì)結(jié)果.本文與該文獻(xiàn)所提出的格式的不同之處在于:
·本文的分析同時(shí)適用于二維與三維情形;
·本文除了適用于單純形網(wǎng)格,還適用于任意多邊形(多面體)網(wǎng)格,文獻(xiàn)[29]只針對二維的三角形網(wǎng)格;
·本文采用Crack-Nicolson格式,時(shí)間方向上精度為二階,文獻(xiàn)[29]對時(shí)間的離散采用了向后Euler差分格式,時(shí)間方向上精度為一階.
本文的剩余內(nèi)容安排如下: 第2節(jié)給出本文的記號和預(yù)備性結(jié)果;第3節(jié)給出半離散格式以及相應(yīng)的誤差估計(jì);第4節(jié)給出了全離散格式以及相應(yīng)的誤差分析;第5節(jié)給出數(shù)值結(jié)果.
引入空間
L(0,T;Hs(D)):={v:(0,T]→Hs(D),v是可測函數(shù),且<}
及相應(yīng)的范數(shù)
‖u‖L(0,T;Hs(D))
在本文中,除特殊說明外,我們用C表示與網(wǎng)格尺度h和時(shí)間步長Δt(稍后定義)無關(guān)的正常數(shù).
Gronwall不等式: 設(shè)u(t),β(t),α(t)是定義在區(qū)間[a,b]上的實(shí)值連續(xù)函數(shù),其中β(t)≥0對任意t∈[a,b],若
則有
?t∈[a,b]
(2)
離散的Gronwall不等式: 設(shè)(kn)n≥0,(pn)n≥0均為非負(fù)數(shù)列.若數(shù)列(φn)n≥0滿足
則
(3)
右矩形公式誤差估計(jì): 設(shè)實(shí)數(shù)a,b滿足a
(4)
梯形公式誤差估計(jì): 設(shè)實(shí)數(shù)a,b滿足a
(5)
則
(6)
根據(jù)文獻(xiàn)[21],我們來定義離散的弱梯度算子.給定Ω上的一個(gè)劃分Th,對于T∈Th,用V(T)表示T上的弱函數(shù)空間,其定義如下:
V(T)={v={v0,vb}:v0∈L2(T),
vb∈H1/2(?T)}
(7)
選擇一個(gè)有限維向量空間G(T)?H(div,T),定義離散弱算子w,G,T:V(T)→G(T).
定義3.1對于任意的v∈V(T),相應(yīng)的離散弱梯度算子w,G,Tv∈G(T) 滿足下列方程:
〈vb,τ·nT〉?T?τ∈G(T)
(8)
特別地,如果G(T)=[Pr(T)]d,我們直接記w,rw,G.
vhb|e∈Pk-1(e),?T∈Th,?e∈εh},
(9)
初值條件為uh(x,0)=Ehu0(定義見后文式(29)),這里
a(uh,vh)=(w,k-1uh,w,k-1vh),
as(uh,vh)=a(uh,vh)+s(uh,vh).
我們引入一種半范數(shù):對于任意vh∈Vh,定義
‖|vh|‖2:=as(uh,uh)=‖wvh‖2+
(10)
引理3.2[21]對于qh={vh0,vhb}∈Vh,T∈Th,wvh=0 在T上成立當(dāng)且僅當(dāng)vh0=vhb=常數(shù).
引理3.3[3]設(shè)整數(shù)m滿足1≤m≤j+1.則對于任意T∈Th,e∈∈h,有
?v∈Hm(T),
?v∈Hm(T),0≤s≤m,
?v∈Hm(T),0≤s+1≤m.
設(shè)uh0(t)|T=Φ0α(t),uhb(t)|T=Φbβ(t),wuh(t)|T=Φγ(t),其中給出一些矩陣的定義如下:
在式(8)中取τ=φj,j=1,2…r3,可得
(11)
(12)
在式(9)中取vh={0,vhb}={0,φjb},j=1,2…r2,可得
(13)
設(shè)
則
(14)
(15)
(16)
(17)
(18)
其中,
在方程(17)和(18)兩端對時(shí)間變量求導(dǎo),則有
設(shè)
0≤t≤T,
(19)
引理3.5[28]設(shè)u∈L(0,T;Hk+1(Ω))為問題(1)的真解. 記ρ{ρ0,ρb}=Ehu-Qhu.則當(dāng)h足夠小時(shí),有
注1文獻(xiàn)[28]中只證明了二維RT元的情形.事實(shí)上,對(k,k-1,k-1)型元以及三維的情況可類似證明.
對(19)兩邊關(guān)于t求導(dǎo),再由類似引理3.5的方式我們可以得到ρt,ρtt的估計(jì)如下.
引理3.6假設(shè)引理3.5的條件均成立,且ut∈L(0,T;Hk+1(Ω)),則當(dāng)h足夠小時(shí),有
進(jìn)一步,若utt∈L(0,T;Hk+1(Ω)),則當(dāng)h足夠小時(shí),有
引理3.7[21]假設(shè)v∈Hk+1(Ω),其中k≥1.則‖v-w(Qhv)‖≤Chk‖v‖k+1,其中
(20)
(21)
(22)
由Ehu的定義可知
(23)
結(jié)合(9),(22)和(23)式可得
((ξ0)t,vh0)+as(ξ,vh)=
(24)
取vh=ξt.由于
則
(25)
對時(shí)間變量從0到t積分,考慮到ξ(x,0)=0,則有
(26)
對于右端第二項(xiàng),由分部積分公式可得
結(jié)合‖|·|‖的定義,有
利用Gronwall不等式,結(jié)合引理3.6,可得如下結(jié)果:
則式(20)可通過引理3.5,引理3.7和三角不等式得到.
在式(24)中取v=ξ,有
C‖(ρ0)t‖2+‖ξ0‖2+
易得
對時(shí)間變量從0到t積分,可得
利用Gronwall不等式,考慮到ξ(x,0)=0,結(jié)合引理3.5,可得如下結(jié)果
結(jié)合引理3.5,引理3.4和三角不等式,我們得出式(21).證畢.
vh)=(f(·,tn+1/2),vh0)
(27)
定理4.1問題(1)的的全離散格式(27)有唯一解.
在開始分析誤差之前,我們先給出一些逼近誤差的定義和估計(jì).當(dāng)n=0,1…N-1時(shí),記
Eh0ut(·,tn+1/2).
當(dāng)n=1,2…,N-1時(shí),記
(28)
(29)
(30)
(31)
根據(jù)Eh0u(·,tn+1)和Eh0u(·,tn)在tn+1/2處的Taylor展開式,可得
由復(fù)合梯形公式(7)和右矩形公式(4)的誤差估計(jì),得
由梯形公式的誤差估計(jì)(5),有
(32)
由引理3.6,引理3.4,3.7和三角不等式,可得
(33)
類似可得
進(jìn)而可得出式(28)~(31).
注2證明過程中利用到了時(shí)間與空間的獨(dú)立性,即(wEhu)tt=wEhutt,(wEhu)t=wEhut,(Eh0u)tt=Eh0utt.
定理4.3設(shè)u是問題(1)的解,且滿足u∈L(0,T;Hk+1(Ω)),utt1,2…N)是全離散格式(26)的解,則h足夠小時(shí),有
(34)
(35)
證明 記
Qhu(x,tn+1/2),n=0,1…,N-1.
(36)
在式(36)中取t=tn+1/2,結(jié)合式(27),(19),可得
(37)
取vh=ξn+1+ξn.我們有
(38)
其中
對于式(38)左端第一項(xiàng),通過計(jì)算可得
(39)
下面估計(jì)式(38)的右端項(xiàng).利用Young不等式,有
(40)
(41)
利用Young不等式,可得
進(jìn)而有
(42)
結(jié)合式(38)~(42)以及|‖·‖|的定義,可得
上式兩邊同時(shí)乘上Δt,對n從0 取到m-1求和,由
再結(jié)合ξ0=0可得
其中
利用離散的Gronwall不等式可得
由引理3.5可得
再結(jié)合引理4.2有
結(jié)合引理3.5,引理3.7與三角不等式,我們得出式(34).
在式(37)中取vh=?tξn+1.對于 0≤n≤N-1,有
(43)
其中
對于方程(43)的左端項(xiàng),有
as(ξn+1+ξn,ξn+1-ξn)=|‖wξn+1‖|2-
(44)
對于其右端項(xiàng),我們有
(45)
同時(shí),由
(46)
結(jié)合式(43)~(45)可得
(47)
上式對n從 0 取到m-1求和,再兩邊同時(shí)乘Δt,可得
(48)
現(xiàn)分別估計(jì)上式右端的三項(xiàng).我們有
對于第三項(xiàng),采用與式(42)同樣方式可得
因此
(49)
結(jié)合式(48),(49)和 |‖·‖|的定義,有
其中
利用離散的 Gronwall 不等式,可得
由引理3.5和引理4.2,我們有
Hk+1(Ω))].
最后,根據(jù)引理3.5,引理3.4和三角不等式,即得式(35).證畢.
本節(jié)將用數(shù)值算例來驗(yàn)證定理4.3中的理論結(jié)果.在所有算例中,取Ω=(0,1)×(0,1),T=1,選擇右端項(xiàng)f使初值條件使真解為u(x,t)=t2x1(1-x1)x2(1-x2).
由表1和表2可以看出,k=1時(shí),u的逼近在時(shí)間精度上有2階,在空間精度上有2階,與定理4.3相符;u的逼近在時(shí)間精度上有2階,在空間精度上有1階,也與定理4.3相符.
表1 正方形網(wǎng)格,k=1Tab.1 Squaredgnid,k=1
表2 正方形網(wǎng)格,k=1Tab.2 Squaredgnid,k=1