一類偏微分方程的差分解法

袁毅楓

【摘要】本文提出把一類橢圓區域上的偏微分方程變換成圓區域上的偏微分方程，然后在極坐標下作網格劃分，使所有網格邊界點與邊界重合，并給出相應的差分方程和解法，從而避免了一般差分解法由邊界條件轉移問題所產生的誤差.

【關鍵詞】差分；方程；解

一、問題的提出

考慮方程：

λ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，….

兩族直線的交點為網格點.

Dk={（xi，yj）|（xi，yj）∈D，并且（xi+1，yj），（xi-1，yj），（xi，yj+1），（xi，yj-1）∈D}為內點.

Dh={（xi，yj）|（xi，yj）∈D-Dh}為網格邊界點.

顯然Dh與D不重合.

在Dh上建立差分格式：

Lhuij=f（xi，yj），（xi，yj）∈Dh，

uij=g（xi，yj），（xi，yj）∈Dh.

（Lh為差分算子）于是出現了邊界條件從D到Dh轉移的問題，帶來了誤差并逐層傳播，影響所得結果的準確性，所以考慮對這類問題做某種變換，同時采取某種特別的解域劃分，使網格邊界點與邊界重合，并給出這種劃分下相應的差分方程和解法.

二、建立差分方程

令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）的系數無定義，因此，只有在r>0的情形下其解才有意義，為了定出有定義的解，需要補充在r=0處u有界.

沿徑向r和角度θ方向取等步長Δr和Δθ作網格劃分：

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.

利用循環邊界條件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個網格點上的差分方程.

在J個網絡邊界點上u（1，θj）=α2（θj），

j=0，1，2，…，J-1.（9）

至此得到了方程（4）的一個差分格式，并且使得網格邊界點與邊界點完全重合.

三、求解差分方程

把（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）

網格邊界點上ui，j由（9）代入.

把D2的內點按如圖所示的排列給予序號.

于是可得方程組（10）的矩陣形式

AU=F，（11）

其中A為方程組的系數矩陣，U為uij依序排列的矩陣，F為f2（ri，θj）依序排列的矩陣.

為了證明所得差分方程是可解的，需要用到下述定理.

極值定理：設差分算子

Lhuij=αijuij-∑blkulk，（xl，yk）∈G（xi，yj），

其中G（xi，yj）表示格式中除（xi，yi）之外相鄰網格點集合，在Dh∪Dh上，算子的系數滿足

aij>0，（xi，yj）∈Dh，

blk>0，（xi，yj）∈G（xi，yj），

dij=aij-∑blk≥0.

uij是定義在Dk上的函數，uij≠常數，如果Lhuij≤0，（xi，yj）∈Dh，則uij不可能在內點取正的最大值：

maxuij≤maxuij.

Dh Dh

如果Lhuij≥0，（xi，yj）∈Dh，則uij不可能在內點取得負的最小值：

minuij≤minuij.

Dh Dh

證明略.

可以驗證差分方程（10）滿足極值定理，由此可得差分格式的穩定性和收斂性.另外，由極值定理可知，滿足差分方程（10）的uij除非在所有網格點上均取同一數值，不會在網格區域的內點取到最大值和最小值.由此立刻可以證明，相應于α2=0情形的差分方程，即在邊界點取零值的差分方程，只有全為零的解，這說明相應于線性方程組（11）的齊次方程組只有零解，因此，線性代數方程組（11）必有唯一解，這說明了所得的差分方程恒有唯一解.

由此可以采用直接的方法來求解，但由于在實際計算中為了保證精確度，步長往往取得很小，網格點很多，所得的差分方程是一個高階線性代數方程，而用直接求解法就會因為受計算機存儲量的限制而產生困難.

由于每個內點上的值只依賴周圍相鄰四點上的值，相應地，線性方程組的系數矩陣有大量零元素，即是一個稀疏陣，采用迭代法求解就可以充分利用這個優點.為此這里給出如下辦法求解這個差分方程.

首先任意給定在網格區域內點（xi，yi）上的數量{u（0）ij}及{u（1）i-1}（i=0，1，2，…，I-1）作為解的近似，把這組數值代入迭代式

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

其中ω是一個常數，可在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為網格點的總數小于某個預先給定的適當小的控制數ε>0時，就可以結束迭代過程，完成差分方程的求解.