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

(廣西師范大學 數(shù)學與統(tǒng)計學院，廣西 桂林 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]等。目前，針對方程中的非線性項，學者們主要通過隱顯格式(IMEX)或線性化方法進行處理，其原因在于全顯式時間離散方法對時空步長要求高，全隱式時間離散方法則在每一個時間點都要迭代求解一個非線性隱式方程組，而隱顯格式或線性化方法則避免了求解非線性代數(shù)方程組，從而減少計算量并提高計算效率。如Shampine等[17]提出一種求解剛性擴散反應偏微分方程的隱顯方法；而后，Hundsdorfer等[18]對具有一般單調(diào)性和有界性的線性多步方法進行了隱顯格式的推廣；Avazzadeh等[1]借助梯形公式使積分線性化表示；Xiao等[19]、Zhang等[20-21]考慮了隱顯單支方法求解幾類時滯微分方程，克服解的時間(偏)導數(shù)具有間斷性所導致的困難，獲得了方法的收斂性結果；Li等[22]、Zhang等[23]則利用泰勒展開等方式對非線性的非齊次項進行了相應的處理。

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

(4)

(5)

(6)

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

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

1 緊隱顯BDF方法格式

(7)

先引入一個常用結果。

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

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

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

(8)

由泰勒展開可知

(9)

(10)

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

(11)

利用引理1可得

(12)

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

(13)

式中

將已知初邊值條件(2)～(3)對應改寫為

(14)

(15)

(16)

定理1假設條件(4)～(6) 成立，則系統(tǒng)(15)～(16)的局部截斷誤差滿足

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。事實上，矩陣A是對角占優(yōu)矩陣，所以是非奇異的，故方程(18)存在唯一解。

3 收斂性分析

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

引理3[23]對任意v∈V，可得

即

(19)

接下來，用數(shù)學歸納法對收斂性定理進行證明。將式(19)中的項拆分討論，對左邊第1項有

(20)

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

(21)

對右邊第1項，利用赫爾德不等式可得

(22)

對于右邊第2項，可以進行如下估計

(23)

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

將上述不等式兩邊同時乘以2τ，并對k求和，可得

再利用引理2的Gronwall不等式，上式可轉化為

最后，應用引理3，得

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

4 數(shù)值算例

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

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

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

在實際驗算過程中，程序的計算運行時間均在3 s內(nèi)，這確保了算法的高效性。

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

此方程的解析解為u(x,t)=e-xt。表1給出了同一空間位置上不同時間點的準確值，而后給出了τ=0.01與τ=0.002 5條件下的截斷誤差?？梢钥闯?，網(wǎng)格越細，所得誤差越小，當空間層取點步長減半時，其誤差減少了1/16，驗證了收斂可達4階。

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

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

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

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

圖3 t=0.5與t=1時近似解與真解的對比

例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,

并在非齊次項中有g(x,s)=ex+tu(x,s), 及

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

此例的解析解為u(x,t)=(1-x6)sin(x+t)。表2詳細給出了在空間點x=0.5處部分時間節(jié)點上的準確值、收斂階以及τ=0.01與τ=0.002 5條件下的截斷誤差?？梢钥闯?，最終的收斂階仍然接近4階。

表2 例2相關的數(shù)值結果

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

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

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

圖6 t=0.5與t=1時近似解與真解的對比

5 結束語

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