陳貞忠 ,馬小霞
(1.新鄉學院數學系,河南新鄉 453000;2.焦作大學基礎部,河南焦作 454003)
近期出現了許多求解拋物型偏微分方程的差分格式[1-8],但交替方向法仍然是求解二維和三維拋物型方程比較理想的數值方法。該法的最早提出者是Peaceman和Rachford,他們提出了二維問題的PR格式[9],此后,又產生了適用于二維和三維問題的 Douglas格式[10],所有這些格式的精度都較低,截斷誤差僅為O(τ2+h2)。
本文對三維常系數非齊次拋物型方程的初邊值問題

導出了一個高精度恒穩定的緊交替方向差分格式,截斷誤差階達到O(τ2+h4)。然后,使用外推算法得到了O(τ3+h6)階精度的近似解,實驗結果與理論分析完全相符。
將求解區域進行網格剖分。設τ=△t為時間步長;h=△x=△y=△z為空間方向步長;為在節點(jh,kh,lh,nτ)處的網函數值,方程(1)的解函數為u(x,y,z,t)。
記u(jh,kh,lh,nτ)=u(j,k,l,n),由Taylor展開可得:

由式(4)和式(5)得(1-▽t)-1=exp(τ),從而 τ=-ln(1-▽t),即:

在(j,k,l,n+1)處考慮方程(1)有:


構造如下4個差分算子:


取η0=,則有:

故當η0=時,式(9)成立。取η1=,則有:

同理取η2=η3=2 ,則有:



式(16)左端可以化為:

式(16)右端可以化為:

由于:

則式(16)可化為:


可得方程(1)的差分格式:

可知格式的截斷誤差為ο(τ2+h4)。
式(19)可分解為:

觀察交替方向算法(20)在每個時間層上只需用追趕法解 3個三對角的方程組,因此計算量較小。
利用Fourier穩定性分析方法,令:

將上式代入格式(20)的誤差方程,即式(20)的齊次形式中,經計算整理,并利用關系式



根據s1的取值范圍可知:

定理 緊差分格式(20)絕對穩定且以ο(τ2+h4)的收斂階收斂。
為了提高格式(20)的數值解的精度,使用Richardson外推法,根據前邊的論述知:


在區域D:{0≤x,y,z≤1,t≥0}上對初邊值問題:

用本文格式(20)和三維問題的Douglas格式[11]求數值解,并與精確解u(x,y,z,t)=e-4tsin(x+y+z)相比較,取h=1/10;τ=rh2=r/100;r=1/2,1,計算到n=200時的結果見表1。

表1 各種算法計算結果與精確解數值比較表
由表1可以看出:本文格式(20)解與精確解均有較好的吻合,它較文獻[11]中的Douglas格式至少精確 2位有效數字,本文格式(20)外推一次所得數值結果與精確解至少有 9位有效數字吻合,這與理論分析一致。
[1] 曾文平.多維拋物型方程的分支絕對穩定的顯格式[J].高等學校計算數學學報,1997,19(2):112-121.
[2] 馬明書.解三維拋物型方程的一個新的高精度顯格式[J].應用數學和力學,1998,19(5):65-469.
[3] 孫鴻烈.解高維熱傳導方程的一族高精度顯式差分格式[J].高校應用數學學報:A輯,1999,14(4):427-432.
[4] 馬明書,王同科.三維拋物型方程的一族高精度分支差分格式[J].應用數學和力學,2000,21(10):1087-1092.
[5] 劉繼軍.二維熱傳導方程的三層顯式差分格式[J].應用數學和力學,2003,24(5):537-543.
[6] 李雪玲,孫志忠.二維變系數反應擴散方程的緊交替方向差分格式[J].高等學校計算數學學報,2006,28(1):83-95.
[7] 馬菊意,楊輝.二維拋物線方程的一個高精度PC格式[J].工程數學學報,2008,25(2):373-376.
[8] 馬明書,王曉峰,馬文娟.二維變系數非齊次拋物型方程的緊交替方向差分格式[J].河南科技大學學報:自然科學版,2008,29(3):88-90.
[9] Peaceman DW,Rachford JH H.The Numerical Solution of Parabolic and Elliptic Differential Equations[J].J Soc Ind Appl Math,1959,3:28-41.
[10] Douglas JJ.Alternating Direction Methods for Three Space Variables[J].Numer Math,1961(4):41-63.
[11] 胡健偉,湯懷民.微分方程數值解[M].北京:科學出版社,1999:229-240.