擴(kuò)散方程Crank-Nicolson格式的穩(wěn)定性

余昌彪,郭 紅,杜明洋,張?zhí)鹛?許文文

齊魯工業(yè)大學(xué)(山東省科學(xué)院) 數(shù)學(xué)與統(tǒng)計學(xué)院,濟(jì)南 250353

關(guān)于擴(kuò)散方程Crank-Nicolson格式的研究已有大量結(jié)果,文獻(xiàn)[1]構(gòu)造了擴(kuò)散系數(shù)為1時的一維和二維拋物方程的Crank-Nicolson格式,并證明它是階數(shù)為二的無條件穩(wěn)定的差分格式;文獻(xiàn)[2]把二維Crank-Nicolson格式,由常系數(shù)推廣到變系數(shù)情形,并證明了它是階數(shù)為二的無條件穩(wěn)定的差分格式。穩(wěn)定性作為數(shù)值格式的一種屬性,在某種意義下,格式的穩(wěn)定性與一致性相結(jié)合,就可以推出數(shù)值格式對于原問題的估計的收斂性。而一致性易于驗證,從而研究穩(wěn)定性是得到收斂性的主要工作。P.D.Lax借助泛函分析的工具,對一類很廣的初值問題證明了差分格式解的收斂性和穩(wěn)定性的等價定理[3]。J.Douglas采用初等的方法討論了具體問題,得到了收斂速度的估計[4]。李榮華綜合了Lax和Douglas工作的優(yōu)點(diǎn),給出Lax等價定理的一個初等證明,去掉了關(guān)于初值問題的適定性要求和簡化了Douglas的證明[5]。而本文研究的任意正常數(shù)的擴(kuò)散系數(shù)的Crank-Nicolson差分格式不僅是無條件穩(wěn)定的,而且誤差較小,計算速度快。

1 擴(kuò)散方程的解析解

考慮常系數(shù)擴(kuò)散方程

(1)

其中a為正常數(shù)且g(x)不恒為零。

令u(x,t)=h(x)p(t)且是方程(1)的充分光滑解,所以u(x,0)=h(x)p(0)=g(x),由于g(x)不恒為零,那么h(x)=g(x)/p(0),進(jìn)而

由擴(kuò)散方程(1)得

它的解是

其中c是任意常數(shù),故擴(kuò)散方程的解析解為

(2)

從u(x,t)的表達(dá)式可以看出擴(kuò)散方程(1)解析解的正負(fù)性完全由g(x)的正負(fù)性控制。

2 差分格式的建立

將區(qū)域[0,1]×(0,+∞)沿x軸和t軸方向進(jìn)行矩形剖分,其中空間步長為Δx=h=1/J,時間步長為Δt=τ,網(wǎng)格點(diǎn)(xj,tn)記作

xj=jΔx=jh,j=0,1,…,J,

tn=nΔt=nτ,n=0,1,2,…。

(3)

初值條件和邊值條件離散為

(4)

(5)

且該格式的截斷誤差為O(τ2+h2),當(dāng)τ和h都趨于零時,O(τ2+h2)也趨于零,所以該差分格式滿足相容性。

|λj(G(τ,k))|≤1+Mτ,j=1,2,…,p,

(6)

其中，λj(G(τ,k))表示增長矩陣G(τ,k)的特征值；M為常數(shù)。條件(6)被稱為von Neumann條件。

引理2:如果差分格式的增長矩陣G(τ,k)是正規(guī)矩陣,則von Neumann條件是差分格式穩(wěn)定的充要條件。

引理3:如果差分格式un+1=Aun的矩陣A是一個正規(guī)矩陣,則譜半徑條件

ρ(A)≤1+Mτ

3 Fourier方法

定理1:若a為正常數(shù)且gj在網(wǎng)格點(diǎn)上有意義,存在常數(shù)τ0>0,使得當(dāng)τ≤τ0,nτ≤T時,則對所有k∈R,擴(kuò)散方程(1)的Crank-Nicolson格式(3)-(5)是無條件穩(wěn)定的。

證明:將(3)式改寫為

(7)

消去公因子有

(2(1+aλ)-aλeikh-aλe-ikh)vn+1=(2(1+aλ)+aλeikh+aλe-ikh)vn,

化簡得

(1+aλ-aλcoskh)vn+1=(1-aλ+aλcoskh)vn,

由此得增長因子

G(τ,k)=(1-aλ+aλcoskh)/(1+aλ-aλcoskh),

所以有

因為a>0,λ>0,所以|G(τ,k)|≤1,那么|G(τ,k)n|≤1,又Crank-Nicolson格式是常系數(shù)差分格式,由引理1和引理2知Crank-Nicolson格式(3)-(5)是無條件穩(wěn)定的。

4 矩陣方法(直接方法)

定理2:若a為正常數(shù)且gj在網(wǎng)格點(diǎn)上有意義,存在常數(shù)τ0>0,使得當(dāng)τ≤τ0,nτ≤T時,則擴(kuò)散方程(1)的Crank-Nicolson格式(3)-(5)是無條件穩(wěn)定的。

證明:可以把(7)式寫為向量形式,即

=

(8)

若令

Aun+1=Bun,

(9)

其中

令C=A-1B,則(9)式可以寫為

un+1=A-1Bun=Cun。

令J-1階方陣

則A和B可以表示為

其中E是J-1階單位矩陣。

又因為矩陣S的特征值[6]為2coskhπ,所以A的特征值為1+aλ-aλcoskhπ,B的特征值為1-aλ+aλcoskhπ,從而矩陣C的特征值為

又矩陣C為對稱矩陣,從而C為正規(guī)矩陣,所以對任意的aλ矩陣C的譜半徑都有

由于問題是線性的,故由引理3知Crank-Nicolson格式(3)-(5)是無條件穩(wěn)定的。

5 數(shù)值例子

為驗證式(3)-(5)的穩(wěn)定性,考慮初邊值問題

(10)

根據(jù)式(2)可知解析解為

u(x,t)=e-π2tsinπx,0≤x≤1,t≥0。

取h=0.02,τ=0.002 5,λ=τ/h2=6.25為網(wǎng)格比。對不同的系數(shù)a,用Crank-Nicolson格式求出擴(kuò)散方程(10)在(0.4,0.02)數(shù)值解,并將它們與解析解的值加以比較。

表1給出了當(dāng)網(wǎng)格比λ=6.25時不同系數(shù)a下的精確解和數(shù)值解,從表中不僅能看出誤差較小,同時也驗證了格式滿足2階的收斂精度。

表1 數(shù)值解與解析解間的比較

從圖1和圖2可以看出,式(3)-(5)的數(shù)值解與精確解的吻合度很好。以上結(jié)論均表明Crank-Nicolson格式是無條件穩(wěn)定的。

圖1 aλ=0.2時的數(shù)值解和解析解

圖2 aλ=2.0時的數(shù)值解和解析解