具有變系數(shù)四階拋物方程的B樣條有限元法

秦丹丹, 唐鑫鑫, 黃文竹

(1. 空軍航空大學(xué) 基礎(chǔ)部, 長春 130022； 2. 吉林大學(xué) 數(shù)學(xué)學(xué)院, 長春 130012；3. 貴州醫(yī)科大學(xué) 生物與工程學(xué)院, 貴陽 550025)

微分方程是一種常見的數(shù)學(xué)模型. 由于實際應(yīng)用中很多問題不僅與空間變量有關(guān), 而且受時間變量影響, 因此拋物方程尤其是高階拋物方程在許多領(lǐng)域應(yīng)用廣泛. 如Cahn-Hilliard方程[1-3]被人們用于分析自然界中的擴散現(xiàn)象[4-5], 其形式如下：

ut+γΔ2u-Δf(u)=0,

其中f(u)在不同的問題中表達式不同.

由于很多微分方程無法找到解析解, 因此, 人們致力于研究微分方程的數(shù)值解法, 其中有限元方法是一種有效的方法. 關(guān)于Cahn-Hilliard方程的有限元法目前已有很多研究成果[6-11]. 針對帶有梯度能量系數(shù)的黏性Cahn-Hilliard方程, Choo等[9]構(gòu)造了有限元格式并進行了理論分析; Qin等[11]討論了晶體表面生長模型的B樣條有限元方法及其近似解的有界性和收斂性, 并通過數(shù)值結(jié)果驗證了格式的有效性. 有限元方法本質(zhì)上是利用分片多項式函數(shù)離散連續(xù)問題, 實現(xiàn)從無限維空間到有限維空間的轉(zhuǎn)化, 因此, 有限維子空間基函數(shù)的選取至關(guān)重要. B樣條作為一種分段多項式函數(shù)[12], 用于各種數(shù)值計算問題, 包括有限元分析[13-15]. 由于B樣條具有緊支集且只有一族基函數(shù), 可節(jié)省存儲空間, 提高計算效率, 因此本文選取B樣條作為有限元基函數(shù).

1 模 型

本文取f(u)=-u, 并限定初邊值條件, 為使模型的適用范圍更廣, 將常系數(shù)γ換成變系數(shù), 得到模型:

(1)

這里I=[0,1]. 變系數(shù)α(x,t)滿足連續(xù)性:

(2)

有界性:

0

(3)

以及關(guān)于時間變量偏導(dǎo)數(shù)的有界性:

(4)

其中s,S,M是正常數(shù).

2 B樣條有限元法

基于邊值條件, 定義函數(shù)空間

(5)

下面給出等距B樣條對應(yīng)卷積的定義.

定義1[11]設(shè)m為正整數(shù),N1(x)為[0,1]上的特征函數(shù), 則函數(shù)

稱為m階B樣條.

本文考慮求解模型(1)的三次B樣條有限元法. 定義剖分Ih: 0=x0

近似解uh(x,t)∈Uh可表示成

(6)

為分析B樣條有限元格式解的有界性與收斂階, 引入下列算子.

a(u-Rhu,vh)=(α(x,t)D2(u-Rhu),D2vh)=0, ?vh∈Uh,

(7)

這里a(u,v)=(α(x,t)D2u,D2v).

由α(x,t)的性質(zhì)(3)和橢圓投影算子的定義(7), 可知其為雙線性的對稱正定算子, 且滿足式(7)的Rhu唯一確定, 并有如下估計[17]:

‖u-Rhu‖+h‖u-Rhu‖1+h2‖u-Rhu‖2≤Ch4‖u‖4.

(8)

下面證明半離散問題解的有界性.

‖uh(t)‖2≤C‖uh(0)‖2, 0≤t≤T,

(9)

其中正常數(shù)C僅與s,M,T有關(guān), 與步長h無關(guān).

證明： 問題(6)是一個常微分方程組, 因此, 它在區(qū)間[0,tn)上存在唯一的局部解. 若可以證明式(9)成立, 再結(jié)合拓展定理, 可得整體解的存在唯一性. 所以只需證明式(9)成立.

在問題(6)中, 取vh=uh, 由式(3)和內(nèi)插不等式知,

易見

(10)

由Gronwall不等式, 可得近似解的L2模有界性：

‖uh(t)‖2≤et/s‖uh(0)‖2≤eT/s‖uh(0)‖2≤C‖uh(0)‖2, 0≤t≤T.

(11)

對式(10)關(guān)于時間變量t積分, 得

(12)

將式(11)代入式(12), 得

(13)

在式(6)中, 再取vh=uh,t, 得

‖uh,t‖2+(α(x,t)D2uh,D2uh,t)-(Duh,Duh,t)=0.

通過簡單求導(dǎo)運算, 可得

由變系數(shù)的性質(zhì), 易得

關(guān)于變量t積分, 得

由ε-不等式可知

由式(11),(13), 可得近似解的H2半模有界性:

‖D2uh(t)‖2≤C‖D2uh(0)‖2, 0≤t≤T.

(14)

又因為

所以‖uh(t)‖2≤C‖uh(0)‖2, 這里正常數(shù)C與s,S,M和uh(0)有關(guān). 證畢.

下面給出半離散問題解的誤差分析.

定理3設(shè)u和uh分別是弱形式(5)和半離散問題(6)的解, 且有u(0)∈H4(I),u,ut∈L2(0,T;H4(I)), 則當(dāng)0≤t≤T時, 近似解的L2模有如下誤差估計：

(15)

其中正常數(shù)C僅與α(x,t)有關(guān), 與步長h無關(guān).

證明: 記θ(t)=Rhu-uh,ρ(t)=u-Rhu, 則u-uh=θ(t)+ρ(t), 從而

‖u-uh‖≤‖θ(t)‖+‖ρ(t)‖.

(16)

推導(dǎo)誤差只需要估計‖θ(t)‖. 聯(lián)立式(5)～(7), 得

(θt,vh)+(α(x,t)D2θ,D2vh)=-(ρt,vh)+(Dθ+Dρ,Dvh).

(17)

在式(17)中取vh=θ, 則有

整理可得

(19)

由Gronwall不等式, 可得

由于

‖θ(0)‖=‖u(0)-uh(0)+Rhu(0)-u(0)‖≤‖u(0)-uh(0)‖+‖ρ(0)‖,

(21)

所以當(dāng)0≤t≤T時, 結(jié)合式(20)和式(21)可知式(15)成立. 證畢.

定理4設(shè)u是弱形式(5)的解,uh是半離散問題(6)的解,u(0)∈H4(I),u,ut∈L2(0,T;H4(I)), 則半離散問題解的誤差估計如下:

(22)

這里正常數(shù)C僅與α(x,t)有關(guān), 與步長h無關(guān).

證明： 在式(17)中選取vh=θt, 則

‖θt‖+(α(x,t)D2θ,D2θt)=-(ρt,θt)+(Dθ+Dρ,Dθt).

通過求導(dǎo)運算, 由ε-等式可得

基于α(x,t)的限制(4)和ε-不等式, 可得

整理可得

(23)

在式(23)兩端關(guān)于變量t積分, 有

再利用α(x,t)的限制條件(3), 可得

由Gronwall引理知

(24)

注意到

‖D2θ(0)‖≤‖D2u(0)-D2uh(0)‖+‖D2Rhu(0)-D2u(0)‖,

(25)

聯(lián)立式(24)和式(25), 有

再由式(16)可得式(22). 證畢.

3 全離散格式

(26)

下面討論全離散格式解的收斂性.

‖u(0)-u0h‖≤Ch4‖u(0)‖4,

(27)

則有如下誤差估計:

(28)

其中C是一個僅依賴于α(x,t)的常數(shù), 與時間步長Δt和空間步長h無關(guān).

(?tθn,vh)+(α(x,tn)D2θn,D2vh)=(rn,vh)+(Dθn+Dρn,Dvh).

(29)

令vh=θn, 則有

整理可得

‖θn‖2-‖θn-1‖2+sΔt‖D2θn‖2≤CΔt(‖θn‖2+‖ρn‖2+‖rn‖2).

關(guān)于n求和, 可得

由投影Rh的性質(zhì)和H?lder不等式知

從而可得rn的估計：

(31)

將式(31)代入式(30), 當(dāng)Δt足夠小, 且CΔt<1-δ(0<δ<1)時, 有

利用式(8)和式(27)可知,

‖θ0‖≤‖u(0)-uh(0)‖+‖u(0)-Rhu(0)‖≤Ch4‖u(0)‖4.

4 數(shù)值模擬

考慮如下問題:

(32)

這里α(x,t)=1+xt, 解析解取為u(x,t)=t2[1-cos(2πx)], 計算可得

g(x,t)=2t-16π3t3sin(2πx)-[16π4t3+(16π4-4π2)t2+2t]cos(2πx).

表1 當(dāng)t=1時, 不同空間步長h對應(yīng)的誤差和收斂階

表2 當(dāng)t=1時, 不同時間步長Δt對應(yīng)的誤差和收斂階

綜上所述, 對于一類四階項帶有變系數(shù)的拋物方程, 本文通過構(gòu)造B樣條有限元法, 引入雙調(diào)和橢圓投影算子, 討論了半離散格式解和全離散格式解的收斂性. 數(shù)值算例很好地驗證了理論結(jié)果. 與Hermite有限元法相比, B樣條有限元法能節(jié)省理更大的存儲空間, 運行速度更快, 而且能達到較好的收斂精度.