正交配置方法求解一維擬線性拋物方程

王佩臣，張志維，樊玉環(huán)，張可為

(黑龍江工程學(xué)院 理學(xué)院，哈爾濱 150050)

考慮下面一維擬線性拋物方程的初邊值問題：

(1)

式中：ψ1(t),ψ2(t),f(u,x,t)和φ(x)為已知函數(shù)。

該方程與其他學(xué)科有著密切的聯(lián)系，它的應(yīng)用直接或間接地源自物理、化學(xué)和生物等學(xué)科領(lǐng)域，具有非常廣泛的實(shí)際應(yīng)用價(jià)值。提出一個(gè)高階數(shù)值方法求解擬線性拋物方程是一個(gè)重要的任務(wù)，當(dāng)前求解該方程主要有有限差分法、有限元、譜方法等，但很少有在時(shí)間和空間都比較穩(wěn)定且有高階精度的方法。文中提出一個(gè)新的方法，組合正交配置方法[1-6]和DIRK[7-8]法來求解一維擬線性拋物方程，該方法在時(shí)間和空間都有很高的精度和穩(wěn)定性。

1 DIRK法

文中簡(jiǎn)介求解一階初始值問題的DIRK法，并給出一個(gè)具體的四階A穩(wěn)定方法。

考慮下面一階初始值問題：

y′(t)=f(t,y),y(t0)=y0,t≥t0.

(2)

式中：y∈Rm,f:R×Rm→Rm，q步DIRK法可表示為

(3)

tn=t0+nh(n=0,1,…),其中，h>0為積分步長(zhǎng)。式(2)、式(3)稱為q級(jí)龍格庫塔法，系數(shù)bi,τi,aij(i,j=1,2,…,q)為實(shí)數(shù)，bi為權(quán)系數(shù)，τi為節(jié)點(diǎn)系數(shù)，A=(aij)q×q為系數(shù)矩陣。這些系數(shù)通常表示為

AτbT

它稱為Butcher’s矩陣。取文獻(xiàn)[7]中Butcher’s矩陣如下

141416-1121449+416073+124115024-19413001449-4160a41a42a43141015372 091-8794112 1362 091+8794112 136141015372 091-8794112 1362 091+8794112 13614

其中

對(duì)于DIRK法，計(jì)算ki(i=1,2,3,4)可以用Newton-Raphson法。Newton-Raphson法計(jì)算過程如下：

2 正交配置方法

首先，用正交配置方法離散方程(1)的空間導(dǎo)數(shù)，得到一個(gè)微分方程組，選取Chebyshev-Gauss-Lobatto節(jié)點(diǎn)作為配置節(jié)點(diǎn)，Legendre多項(xiàng)式作為投影多項(xiàng)式。

設(shè)式(1)的解為

(4)

取N個(gè)配置點(diǎn)a=x1

(5)

由式(5)求解αi(t)得

(6)

由式(4)得u(x,t)的一階和二階偏導(dǎo)為

(7)

(8)

把式(6)代入到式(7)、式(8)中得

(9)

(10)

(11)

(12)

下面選用正交函數(shù)作為實(shí)驗(yàn)函數(shù),首先考慮下面[-1,1]區(qū)間的Legendre多項(xiàng)式。

Legendre多項(xiàng)式：

滿足下面性質(zhì)：

1)遞推關(guān)系。

2)Legendre多項(xiàng)式的正交性。

其中

下面做y∈[-1,1]→x∈[a,b]的映射。

令[a,b]區(qū)間的Chebyshev-Gauss-Lobatto節(jié)點(diǎn)

類似上面的討論，令同樣設(shè)gi(x)=Pi-1(x)

(13)

式(13)寫成多項(xiàng)式展開的形式

(14)

由式(14)得在xj點(diǎn)u(x,t)和u(x,t)的一階和二階偏導(dǎo)為

(15)

(16)

(17)

把式(14) 、式(15)、 式(16)寫成矩陣的形式為

(18)

其中

U(t)=(u(x1,t),u(x2,t),…,u(xN,t))T,

B(t)=(β1(t),β2(t),…,βN(t))T,

C=(cij)N×N,D=(dij)N×N,

則由式(18)得Β(t)=Q-1U(t)，

(19)

(20)

由式(20) 化簡(jiǎn)式(1)得

i=2,3,…,N-1.

(21)

其中,u(x0,t)=ψ1(t),u(xN,t)=ψ2(t)，式(21)化為

miNψ2(t)+f(ui,xi,t),i=2,3,…,N-1.

(22)

式(22)滿足下面初始條件

u(xk,t0)=φ(xk),k=1,2，…,N-1.

(23)

用第二部分的DIRK法求常微分方程組式(22)、式(23)，進(jìn)而得到原方程的解。

3 數(shù)值實(shí)例

3.1 測(cè)試實(shí)例[9]

初始條件為

u(x,0)=sin (πx), 0

u(0,t)=u(1,t)=0,t≥0.

該問題的精確解為

u(x,t)=e-tsin (πx).

圖1為精確解曲面, 圖2為在h=0.01,N=11時(shí)的數(shù)值解曲面，圖3為在h=0.01,N=11時(shí)的誤差曲面。

圖1 精確解曲面

圖2 在h=0.01,N=11時(shí)的數(shù)值解曲面

圖3 在h=0.01,N=11時(shí)的誤差曲面

3.2 測(cè)試實(shí)例

初始條件為

u(x,0)=sinh(πx), 0

u(0,t)=0,u(1,t)=e-tsinh(π),t≥0.

該問題的精確解為

u(x,t)=e-tsinh(πx).

圖4為精確解曲面, 圖5為在h=0.01,N=11時(shí)的數(shù)值解曲面，圖6為在h=0.01,N=11時(shí)的誤差曲面。

圖4 精確解曲面

圖5 在h=0.01,N=11時(shí)的數(shù)值解曲面

圖6 在h=0.01,N=11時(shí)的誤差曲面

4 結(jié)束語

組合使用正交配置法和對(duì)角隱式龍格庫塔法求解一維熱傳導(dǎo)方程，使用正交配置法離散空間導(dǎo)數(shù)，用對(duì)角隱式龍格庫塔法求線性解常微分方程組，通過數(shù)值實(shí)例，可以說明該方法有很高的精度和穩(wěn)定性，數(shù)值實(shí)例證實(shí)該方法的有效性。