非線性偏積分微分方程緊隱顯BDF方法的誤差分析

蘭海峰，肖飛雁，張根根，朱 瑞

(廣西師范大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，廣西 桂林 541006)

本文主要研究一類非線性拋物型偏積分微分方程的數(shù)值解法，其方程的基本形式為

(1)

式中(x,t)∈[a,b]×[0,T],初邊界條件定義為

u(x,0)=φ(x),a≤x≤b,

(2)

u(a,t)=ua(t),u(b,t)=ub(t), 0≤t≤T。

(3)

偏積分微分方程的數(shù)值解法已有許多研究，如有限差分方法[1-6]、譜方法[7-11]以及有限元方法[12-16]等。目前，針對(duì)方程中的非線性項(xiàng)，學(xué)者們主要通過(guò)隱顯格式(IMEX)或線性化方法進(jìn)行處理，其原因在于全顯式時(shí)間離散方法對(duì)時(shí)空步長(zhǎng)要求高，全隱式時(shí)間離散方法則在每一個(gè)時(shí)間點(diǎn)都要迭代求解一個(gè)非線性隱式方程組，而隱顯格式或線性化方法則避免了求解非線性代數(shù)方程組，從而減少計(jì)算量并提高計(jì)算效率。如Shampine等[17]提出一種求解剛性擴(kuò)散反應(yīng)偏微分方程的隱顯方法；而后，Hundsdorfer等[18]對(duì)具有一般單調(diào)性和有界性的線性多步方法進(jìn)行了隱顯格式的推廣；Avazzadeh等[1]借助梯形公式使積分線性化表示；Xiao等[19]、Zhang等[20-21]考慮了隱顯單支方法求解幾類時(shí)滯微分方程，克服解的時(shí)間(偏)導(dǎo)數(shù)具有間斷性所導(dǎo)致的困難，獲得了方法的收斂性結(jié)果；Li等[22]、Zhang等[23]則利用泰勒展開等方式對(duì)非線性的非齊次項(xiàng)進(jìn)行了相應(yīng)的處理。

本文在空間方向采用緊致差分方法，時(shí)間方向采用隱顯BDF方法對(duì)方程(1)～(3)進(jìn)行離散，積分項(xiàng)則利用梯形數(shù)值求積公式求解，最終構(gòu)造出一種高效且穩(wěn)定的數(shù)值方法，稱為緊隱顯BDF方法。本方法只依賴所求節(jié)點(diǎn)的臨近點(diǎn)，故最終可得系數(shù)矩陣為三對(duì)角陣的線性代數(shù)系統(tǒng)，并使用Thomas算法求得最終結(jié)果。為此，設(shè)方程的解析解u(x,t)關(guān)于方程(1)～(3)都是充分光滑的，即u(x,t)∈C6,4([a,b]×[0,T])，設(shè)fμ和fν分別是定義在0-鄰域上對(duì)函數(shù)f第1和第2個(gè)元素的一階偏導(dǎo)，0>0?？梢栽O(shè)

(4)

(5)

(6)

式中：gγ表示g函數(shù)中關(guān)于第3個(gè)元素的一階導(dǎo)數(shù)；π1={x∈(a,b),0

本文安排如下：第1章得出求解方程(1)～(3)的緊隱顯BDF方法及相應(yīng)的局部截?cái)嗾`差；第2章給出相應(yīng)數(shù)值格式的矩陣形式，驗(yàn)證其三對(duì)角性，并分析了可解性；第3章分析緊隱顯BDF方法的收斂性；最后，在第4章給出數(shù)值算例以驗(yàn)證數(shù)值格式的準(zhǔn)確性和有效性。

1 緊隱顯BDF方法格式

(7)

先引入一個(gè)常用結(jié)果。

引理1[23]若存在p(x)∈C6[xi-1,xi+1]，則

式中ω∈(xi-1,xi+1)。

在式(7)的基礎(chǔ)上，考慮在節(jié)點(diǎn)(xi,tk+1)處的離散方程，可得

(8)

由泰勒展開可知

(9)

(10)

對(duì)式(8)兩邊同時(shí)使用算子A，可得到緊差分的數(shù)值格式

(11)

利用引理1可得

(12)

再將式(9)、(10)、(12)代入式(11)中，得

(13)

式中

將已知初邊值條件(2)～(3)對(duì)應(yīng)改寫為

(14)

(15)

(16)

定理1假設(shè)條件(4)～(6) 成立，則系統(tǒng)(15)～(16)的局部截?cái)嗾`差滿足

2 緊隱顯BDF方法的矩陣形式

(17)

式中：1≤i≤M-1；1≤k≤N-1。

(18)

其中系數(shù)矩陣定義如下

A，B都是M-1階矩陣，定義在RM-1的列向量有U=(u1,u2,…,uM-1)T，以及

式中k≥1。事實(shí)上，矩陣A是對(duì)角占優(yōu)矩陣，所以是非奇異的，故方程(18)存在唯一解。

3 收斂性分析

設(shè)定義在Ωh上的網(wǎng)格函數(shù)空間V={v|v=(v0,v1,…,vM),v0=vM=0}，對(duì)任意v,ω∈V，其內(nèi)積與范數(shù)定義如下

引理3[23]對(duì)任意v∈V，可得

即

(19)

接下來(lái)，用數(shù)學(xué)歸納法對(duì)收斂性定理進(jìn)行證明。將式(19)中的項(xiàng)拆分討論，對(duì)左邊第1項(xiàng)有

(20)

利用離散的格林公式，可以將式(19)中的左邊第2項(xiàng)進(jìn)行估計(jì)，

(21)

對(duì)右邊第1項(xiàng)，利用赫爾德不等式可得

(22)

對(duì)于右邊第2項(xiàng)，可以進(jìn)行如下估計(jì)

(23)

將式(20)～(23)帶入式(19)，得

將上述不等式兩邊同時(shí)乘以2τ，并對(duì)k求和，可得

再利用引理2的Gronwall不等式，上式可轉(zhuǎn)化為

最后，應(yīng)用引理3，得

由歸納法原理，定理得證。

4 數(shù)值算例

本章將給出數(shù)值算例以驗(yàn)證算法的準(zhǔn)確性和有效性。由式(17)可知，計(jì)算uk+1的值還需用到uk以及uk-1的值。事實(shí)上，只需提供k=1層，即時(shí)間第1層的值即可。為保證空間上仍為2階收斂，利用泰勒展開，設(shè)存在ηi∈(t0,t1)，使得

則時(shí)間上第1層的節(jié)點(diǎn)信息可表示為

式中φ″(xi)是函數(shù)u對(duì)x的2階導(dǎo)函數(shù)并帶入了t0的值，t0=0使得f中第2項(xiàng)的值為0。

在實(shí)際驗(yàn)算過(guò)程中，程序的計(jì)算運(yùn)行時(shí)間均在3 s內(nèi)，這確保了算法的高效性。

例1根據(jù)方程(1)～(3)的設(shè)定，考慮如下方程

此方程的解析解為u(x,t)=e-xt。表1給出了同一空間位置上不同時(shí)間點(diǎn)的準(zhǔn)確值，而后給出了τ=0.01與τ=0.002 5條件下的截?cái)嗾`差。可以看出，網(wǎng)格越細(xì)，所得誤差越小，當(dāng)空間層取點(diǎn)步長(zhǎng)減半時(shí)，其誤差減少了1/16，驗(yàn)證了收斂可達(dá)4階。

表1 例1的數(shù)值結(jié)果

當(dāng)m與n的分割為10×100時(shí)，例1的近似解的圖像如圖1所示；圖2給出了同等分割下近似解與真解的誤差情況；圖3展示了t=0.5與t=1時(shí)近似解與真解的情況。

圖1 例1所求近似解的值(m=10,n=100)

圖2 例1所求近似解與真解的誤差(m=10,n=100)

圖3 t=0.5與t=1時(shí)近似解與真解的對(duì)比

例2考慮式(1)～(3)具有初始條件

u(x,0)=(1-x6)sinx,0≤x≤1,

與邊界條件

u(0,t)=sint,0≤t≤1;u(1,t)=0, 0≤t≤1,

并在非齊次項(xiàng)中有g(shù)(x,s)=ex+tu(x,s), 及

的拋物型偏積分微分方程。

此例的解析解為u(x,t)=(1-x6)sin(x+t)。表2詳細(xì)給出了在空間點(diǎn)x=0.5處部分時(shí)間節(jié)點(diǎn)上的準(zhǔn)確值、收斂階以及τ=0.01與τ=0.002 5條件下的截?cái)嗾`差?？梢钥闯觯罱K的收斂階仍然接近4階。

表2 例2相關(guān)的數(shù)值結(jié)果

圖4展示了例2在m與n的分割為10×100的條件下所得的近似解的圖像；圖5給出了同等密度分割下近似解與真解的誤差；圖6給出了t=0.5與t=1時(shí)的近似解與真解對(duì)比情況。

圖4 例2所求近似解的值(m=10,n=100)

圖5 例2所求近似解與真解的誤差(m=10,n=100)

圖6 t=0.5與t=1時(shí)近似解與真解的對(duì)比

5 結(jié)束語(yǔ)

本文給出一種求解非線性拋物型偏積分微分方程的高階數(shù)值格式，此格式在空間上采用4階緊差分方法，在時(shí)間上基于2階隱顯BDF方法進(jìn)行離散，并利用梯形求積公式處理了非齊次項(xiàng)中的積分項(xiàng)；而后證明了此方法的收斂性；最后，數(shù)值實(shí)驗(yàn)的結(jié)果表明此方法具有良好的精度，收斂階能夠達(dá)到預(yù)估，驗(yàn)證了方法的高效性和穩(wěn)定性。今后考慮結(jié)合交替方向等技巧將此結(jié)果推廣至2維和3維情形。