無網(wǎng)格法求解一類分段連續(xù)型延遲偏微分方程

鐘 霖， 馬淑芳， 萊 蒙

(東北林業(yè)大學 理學院，哈爾濱 150040)

1 引 言

Myshkis[1]最先發(fā)現(xiàn)存在分段連續(xù)型的延遲微分方程，Wiener等[2,3]對該類方程解的性質進行了研究。該類方程同時具有微分方程和差分方程的性質，能更精確地描述反饋控制和混沌系統(tǒng)等領域的一些問題，因此受到學者的高度重視。延遲微分方程的解析解不易獲得，因此，發(fā)展該類方程的數(shù)值解研究十分必要。目前，對于該類方程數(shù)值解的穩(wěn)定性、振動性和周期性已有大量研究。文獻[4]研究了生態(tài)學中自變量分段連續(xù)型延遲Logistictic方程的數(shù)值解穩(wěn)定性。文獻[5]討論了非線性分段連續(xù)型延遲偏微分方程數(shù)值解的漸近穩(wěn)定性，給出了non-confluent Runge -kutta方法針對標量方程的穩(wěn)定充要條件。文獻[6]運用θ-方法和單腿θ-方法討論了分段連續(xù)型延遲微分方程數(shù)值解的穩(wěn)定性。

近年來，分段連續(xù)型延遲偏微分方程也引起了廣泛的關注。文獻[7]運用Galerkin法研究了分段連續(xù)型延遲偏微分方程，證明了在解析解漸近穩(wěn)定的條件下，半離散和全離散方法的數(shù)值解都是漸近穩(wěn)定的。文獻[8]利用θ-方法研究了超前滯后型分段連續(xù)型偏微分方程,得到了數(shù)值解漸近穩(wěn)定的條件。

本文考慮如下分段連續(xù)型偏微分方程：

(1)

式中t>0,a,b∈R,u(x,t):Ω=[0,1]×[0,∞)→R，v:[0,1]→R，[·]表示最大取整函數(shù)。方程(1)為對流-擴散方程，是一類基本的運動方程。該類方程可以描述在以速度為b運動的河流中污染物濃度u(x,t)的變化,其中a2ux x為擴散項，bux x(x,[t])為離散時間的對流項[11]。文獻[9]采用θ-格式求解方程(1),給出了格式依賴于網(wǎng)格比的穩(wěn)定性條件。當θ=1/2(Crank-Nicolson格式)，網(wǎng)格比是500時，方法是不穩(wěn)定的。這與經(jīng)典的Crank-Nicolson 格式具有無條件穩(wěn)定性不一致。為了克服網(wǎng)格比對數(shù)值方法穩(wěn)定性的影響，本文考慮采用無網(wǎng)格法求解方程(1)。無網(wǎng)格法是建立差分格式時，不需利用預定義的網(wǎng)格信息進行離散的方法[10]，可方便地模擬復雜形狀的流場，并在處理網(wǎng)格移動和畸變等問題時有明顯的優(yōu)勢。

本文主要目的是利用無網(wǎng)格法求解方程(1)，并分析方法的穩(wěn)定條件，最終得到不依賴于網(wǎng)格比的穩(wěn)定性條件。

2 分段連續(xù)型延遲偏微分方程的求解

為得到該方程的離散格式，令時間方向的步長為Δt=1/m，m為正整數(shù)。時間節(jié)點tn=nΔt(n=0,1,…)。

利用θ-加權有限差分(0≤θ≤1)離散方程(1)有

θ[a2ux x(x,t)+bux x[t]]n + 1

(2)

式中un=u(x,tn)。整理有

un + 1=un+(1-θ)Δt[a2ux x(x,tn)+bux x(x,[tn])]+

θΔt[a2ux x(x,tn + 1)+bux x(x,[tn + 1])]

(3)

(4)

即

(5)

在區(qū)間[0,1]上選擇節(jié)點xi(i=1,…,N)。i=2,…,N-1時為內(nèi)點，x1和xN為邊界點。利用徑向基函數(shù)的組合近似u(x,tn)有

(6)

表1 典型的徑向基函數(shù)

將式(6)代入式(5)并配置每個節(jié)點xi(i=2,…,N-1)有

(7)

將以上方程進行簡單的處理有

(8)

通過在邊界點x1和xN處使用邊界條件可得

(9)

(10)

從而，由式(8～10)聯(lián)立可得

(A-θΔta2B)Λk m + l + 1=[C+(1-θ)Δta2B)Λk m + l+

ΔtbBΛk m

(11)

其中

整理可得

DΛk m + l + 1=EΛk m + l+FΛk m

(12)

3 穩(wěn)定性分析

(13)

化簡可得

(l=0,1,2,…,m-1)(14)

(15)

遞歸可得

(16)

整理有

(17)

定理1若|b|≤a2成立。當滿足如下條件時，即

當m為偶數(shù)時

(18)

(19)

方程(11)是穩(wěn)定的。

推論1當θ=1時,該方法是無條件漸近穩(wěn)定的。

注1 式(18,19)中β是傅里葉分析的波數(shù)，其計算公式是β=2π/波長。適當?shù)牟ㄩL可導致式(18,19)的Δt的上界很大，而通常Δt取值很小，故β對于Δt的取值無太大影響。

4 數(shù)值算例

為驗證所得結果，首先給出兩個范數(shù)，

(20)

(21)

算例1[9]考慮方程(22),

(22)

由于本文要研究網(wǎng)格比r=Δt/Δx2=500，θ=1/2時的穩(wěn)定性，故取時間步長為Δt=1/500，取Δx=1/500對應到文獻[9]的網(wǎng)格比為500。由文獻[16]可知，方程(1)在t為整數(shù)時的精確解為[e-a2π2+b(e-a2π2-1)/a2]tsin πx。表2給出了在t=1時，不同c值下，數(shù)值解與精確解的誤差，由表可知總體趨勢是c越小誤差越小，故以下均取c=Δx2。方程(22)的精確解與數(shù)值解誤差列入表3，可以看出，該數(shù)值方法具有較高精度。

表2 不同參數(shù)c下精確解與數(shù)值解的誤差

圖1利用Matlab模擬該方法在t∈[0,10]，Δt=1/500，Δx=1/500，θ=1/2，c=Δx2時的數(shù)值解。可以看出數(shù)值解是漸近穩(wěn)定的。因此解決了文獻[9]網(wǎng)格比為500不穩(wěn)定的現(xiàn)象，證明了該方法的適用性。

算例2[9]考慮方程(23)，

(23)

圖1 方程(22)的數(shù)值解

圖2和圖3為在t∈[0,10]，Δt=1/3000，Δx=1/10，θ=0和0.25時方程(23)的數(shù)值解；圖4和圖5為t∈[0,10]，Δt=1/100，Δx=1/40，θ=0.75和1時方程(23)的數(shù)值解?？梢钥闯觯瑪?shù)值解是漸近穩(wěn)定的,從而驗證了定理1所得的結論。

由于徑向基函數(shù)便于向高維問題推廣，故考慮n維分段連續(xù)型偏微分方程：

(24)

式中U(x,t)∈Rn,V(x)∈R,C0∈Rn,t>0,L,M∈Rn × n。

圖2 方程(23)θ=0時數(shù)值解

圖3 方程(23)θ=0.25時數(shù)值解

圖4方程(23)θ=0.75時數(shù)值解

Fig.4 Numerical solution forθ=0.75 in Eq.(23)

圖5 方程(23)θ=1時數(shù)值解

圖6 數(shù)值解分量圖u 1

圖7 數(shù)值解分量圖u 2

5 結 論

本文采用了MQ徑向基函數(shù)的無網(wǎng)格法求解了分段連續(xù)型延遲偏微分方程，使用了θ-加權有限差分法，然后利用徑向基函數(shù)的組合來近似求解方程，并利用配點法得到了計算數(shù)值解的方程組，最后給出了該方程的穩(wěn)定性分析。數(shù)值算例驗證了方法的有效性和穩(wěn)定性，進一步表明所采用的數(shù)值方法可以有效推廣到n維分段連續(xù)型偏微分方程。