袁毅楓
【摘要】本文提出把一類橢圓區(qū)域上的偏微分方程變換成圓區(qū)域上的偏微分方程,然后在極坐標下作網(wǎng)格劃分,使所有網(wǎng)格邊界點與邊界重合,并給出相應(yīng)的差分方程和解法,從而避免了一般差分解法由邊界條件轉(zhuǎn)移問題所產(chǎn)生的誤差.
【關(guān)鍵詞】差分;方程;解
一、問題的提出
考慮方程:
λ22ux2+2uy2=f(x,y),(x,y)∈D,
u|D=α(x,y),
其中λ≠0,D=(x,y)x2y2+y2<1.(1)
對此類方程的一般差分解法是首先分別沿x,y軸方向取步長為h和k,作兩族與坐標平行的直線:
xi=ih,i=0,±1,±2,….
yj=jk,j=0,±1,±2,….
兩族直線的交點為網(wǎng)格點.
Dk={(xi,yj)|(xi,yj)∈D,并且(xi+1,yj),(xi-1,yj),(xi,yj+1),(xi,yj-1)∈D}為內(nèi)點.
Dh={(xi,yj)|(xi,yj)∈D-Dh}為網(wǎng)格邊界點.
顯然Dh與D不重合.
在Dh上建立差分格式:
Lhuij=f(xi,yj),(xi,yj)∈Dh,
uij=g(xi,yj),(xi,yj)∈Dh.
(Lh為差分算子)于是出現(xiàn)了邊界條件從D到Dh轉(zhuǎn)移的問題,帶來了誤差并逐層傳播,影響所得結(jié)果的準確性,所以考慮對這類問題做某種變換,同時采取某種特別的解域劃分,使網(wǎng)格邊界點與邊界重合,并給出這種劃分下相應(yīng)的差分方程和解法.
二、建立差分方程
令x=λt,(1)成為
2ut2+2uy2=f1(t,y),(t,y)∈D1,
u|D1=α1(t,y),
其中D1={(t,y)|t2+y2<1}.(2)
令t=rcosθ,y=rsinθ,
即r=t2+y2,θ=arctanyt,(2)成為
2ur2cos2θ+ur·sin2θr-2uθr·cosθ·sinθr+2ur2·sin2θ+2uθr·sinθ·cosθr+ur·cos2θr-2urθ·sinθ·cosθr+uθ·sinθ·cosθr2+2uθ2·sin2θr2+2urθ·sinθ·cosθr+2urθ·sinθ·cosθr-uθ·sinθ·cosθr2+2uθ2·cos2θr2-uθ·sinθ·cosθr2=1r·rrur+1r2·2uθ2=f(r,θ),(r,θ)∈D2,
u|D2=α2(1,θ)=α2(θ).
其中D2={(r,θ)|0≤r<1,0≤θ≤2π}.(3)
當r=0時,方程(3)的系數(shù)無定義,因此,只有在r>0的情形下其解才有意義,為了定出有定義的解,需要補充在r=0處u有界.
沿徑向r和角度θ方向取等步長Δr和Δθ作網(wǎng)格劃分:
ri=i+12Δr,i=0,1,2,…,I.
θj=jΔθ,j=0,1,2,…,J-1.
Δr=1I+12,Δθ=2πJ.
于是方程(3)的差分方程為
1ri(Δr)2δr(rδruij)+1r2i(Δθ)2δ2θuij=f2(ri,θj),(4)
其中,δ和δθ分別表示r方向和θ方向的中心差分算子.
即,δruij=ui+12,j-ui-12,j;
δθuij=ui,j+12-ui,j-12.
利用循環(huán)邊界條件u(r,θ)=u(r,θ+2π),具體寫出方程(4)在j=0和j=J-1上的差分方程:
j=0,1ri(Δr2)[ri+12ui+1,0-(ri+12+ri-12)ui,0+ri-12ui-1,0]+1r2i(Δθ)2(ui,1-2ui,0+ui,j-1)=f2(ri,0),
i=1,2,…,I-1.(5)
j=J-1,1ri(Δr)2[ri+12ui+1,j-1-(ri+12+ri-12)ui,j-1+ri-12ui-1,j-1]+1r2i(Δθ)2(ui,0-2ui,j-1+ui,j-2)=f2(ri,θj-1),
i=1,2,…,I-1.(6)
當j=1,2,…,J-2時,差分方程為
1ri(Δr)2[ri+12ui+1,j-(ri+12+ri-12)ui,j+ri-12ui-1,j]+1r2i(Δθ)2(ui,j+1-2ui,j+ui,j-1)=f2(ri,θj),
i=1,2,…,I-1.(7)
令(7)中ri=Δr2,得(r0,θj)處的差分方程:
2(Δr)2(u1,j-u0,j)+4(Δr)2(Δθ)2(u0,j+1-2u0,j-u0,j-1)=f2(r0,θj),
j=0,1,2,…,J-1.(8)
至此建立了I×J個網(wǎng)格點上的差分方程.
在J個網(wǎng)絡(luò)邊界點上u(1,θj)=α2(θj),
j=0,1,2,…,J-1.(9)
至此得到了方程(4)的一個差分格式,并且使得網(wǎng)格邊界點與邊界點完全重合.
三、求解差分方程
把(5)(6)(7)(8)分別整理,并對每處(r0,θj)記為
αi,jui,j-αi+1,jui+1,j-αi-1,jui-1,j-αi,j+1ui,j+1-αi,j-1ui,j-1=-f2(ri,θj),
i=0,1,2,…,I-1;j=0,1,2,…,J-1.(10)
網(wǎng)格邊界點上ui,j由(9)代入.
把D2的內(nèi)點按如圖所示的排列給予序號.
于是可得方程組(10)的矩陣形式
AU=F,(11)
其中A為方程組的系數(shù)矩陣,U為uij依序排列的矩陣,F(xiàn)為f2(ri,θj)依序排列的矩陣.
為了證明所得差分方程是可解的,需要用到下述定理.
極值定理:設(shè)差分算子
Lhuij=αijuij-∑blkulk,(xl,yk)∈G(xi,yj),
其中G(xi,yj)表示格式中除(xi,yi)之外相鄰網(wǎng)格點集合,在Dh∪Dh上,算子的系數(shù)滿足
aij>0,(xi,yj)∈Dh,
blk>0,(xi,yj)∈G(xi,yj),
dij=aij-∑blk≥0.
uij是定義在Dk上的函數(shù),uij≠常數(shù),如果Lhuij≤0,(xi,yj)∈Dh,則uij不可能在內(nèi)點取正的最大值:
maxuij≤maxuij.
Dh Dh
如果Lhuij≥0,(xi,yj)∈Dh,則uij不可能在內(nèi)點取得負的最小值:
minuij≤minuij.
Dh Dh
證明略.
可以驗證差分方程(10)滿足極值定理,由此可得差分格式的穩(wěn)定性和收斂性.另外,由極值定理可知,滿足差分方程(10)的uij除非在所有網(wǎng)格點上均取同一數(shù)值,不會在網(wǎng)格區(qū)域的內(nèi)點取到最大值和最小值.由此立刻可以證明,相應(yīng)于α2=0情形的差分方程,即在邊界點取零值的差分方程,只有全為零的解,這說明相應(yīng)于線性方程組(11)的齊次方程組只有零解,因此,線性代數(shù)方程組(11)必有唯一解,這說明了所得的差分方程恒有唯一解.
由此可以采用直接的方法來求解,但由于在實際計算中為了保證精確度,步長往往取得很小,網(wǎng)格點很多,所得的差分方程是一個高階線性代數(shù)方程,而用直接求解法就會因為受計算機存儲量的限制而產(chǎn)生困難.
由于每個內(nèi)點上的值只依賴周圍相鄰四點上的值,相應(yīng)地,線性方程組的系數(shù)矩陣有大量零元素,即是一個稀疏陣,采用迭代法求解就可以充分利用這個優(yōu)點.為此這里給出如下辦法求解這個差分方程.
首先任意給定在網(wǎng)格區(qū)域內(nèi)點(xi,yi)上的數(shù)量{u(0)ij}及{u(1)i-1}(i=0,1,2,…,I-1)作為解的近似,把這組數(shù)值代入迭代式
u(k+1)ij=ωαi+1αiju(k)i+1,j+αi-1αiju(k+1)i-1,j+αi,j+1αiju(k)i,j+1
+αi,j-1αiju(k+1)i,j-1+(1-ω)u(k)ij,(xi,yi)∈Dk,
u(k+1)ij=α2(θ1),(xi,yi)∈Dk
其中ω是一個常數(shù),可在1<ω<2中選定.當k相當大時,{u(k)ij}就給出所要求的近似解.通常對充分大的k,當相鄰兩次迭代解{u(k-1)ij},{u(k)ij}間的誤差maxi,j|u(k)ij-u(k-1)ij|或1N∑i,j|u(k)ij-u(k-1)ij|,其中N為網(wǎng)格點的總數(shù)小于某個預先給定的適當小的控制數(shù)ε>0時,就可以結(jié)束迭代過程,完成差分方程的求解.