開依沙爾·熱合曼,努爾買買提·黑力力

(新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,新疆 烏魯木齊 830046)

求解擴散方程的二級四階隱式Runge-Kutta方法

開依沙爾·熱合曼,努爾買買提·黑力力

(新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,新疆 烏魯木齊 830046)

對空間變量應(yīng)用中心差分格式和緊致差分格式離散,時間變量采用二級四階Runge-Kutta方法,構(gòu)造求解擴散方程的精度為O(τ4+h2)和O(τ4+h4)的兩種絕對穩(wěn)定的隱式差分格式,討論穩(wěn)定性,并將數(shù)值試驗結(jié)果與Crank-Nicholson格式進行比較,數(shù)值結(jié)果表明該方法是求解擴散方程的有效數(shù)值計算方法之一.

擴散方程;緊致格式;二級四階Runge-Kutta方法;兩層隱格式;Crank-Nicolson格式

0 引言

考慮如下齊次邊界條件的一維擴散方程:

(1)

擴散方程是一類描述物理量隨時間的擴散和衰減規(guī)律的拋物型微分方程.自然環(huán)境、工程設(shè)備及生物機體中的許多物理現(xiàn)象,諸如氣體的擴散、液體的滲透、熱的傳導(dǎo)以及半導(dǎo)體材料中雜質(zhì)的擴散等都可用擴散方程來描述.由于物理問題本身的復(fù)雜性,其精確解往往不容易求得,研究其數(shù)值求解方法具有非常重要的理論意義和工程應(yīng)用價值.

目前求解該問題的差分格式主要有顯式格式、隱式格式、Crank-Nicolson格式[1-3]和其他一些格式[4-11].文獻[4-5]中空間變量用中心差分格式離散化,時間變量分別用四階經(jīng)典R-K方法和修正梯形公式得到了精度分別為O(τ4+h2)和O(τ3+h2)的數(shù)值格式.文獻[6-7]中空間變量用四階緊致差分格式離散化,時間方向用Crank-Nicolson格式構(gòu)造非齊次和齊次擴散方程的絕對穩(wěn)定的差分格式,文獻[8-9]中中空間變量用中心差分格式離散化,時間變量用修正Crank-Nicolson格式得到一維和二維熱傳導(dǎo)方程絕對穩(wěn)定的顯式差分格式.文獻[10]中空間變量用中心差分格式離散化,時間方向用三次C1樣條方法構(gòu)造精度為O(τ4+h2)的絕對穩(wěn)定的兩層隱式格式.文獻[11]在文獻[10]的基礎(chǔ)上對空間變量用四階緊致差分格式離散化,對時間方向用三次C1樣條方法構(gòu)造熱傳導(dǎo)方程和對流擴散方程的精度為O(τ4+h4)的絕對穩(wěn)定的兩層隱式格式.

受上述文獻的啟發(fā),對空間變量分別用中心差分格式和緊致差分格式離散化,時間變量應(yīng)用二級四階隱式Runge-Kutta方法構(gòu)造擴散方程的精度分別為O(τ4+h2)和O(τ4+h4)的兩種絕對穩(wěn)定的隱式差分格式,并討論其穩(wěn)定性,數(shù)值值結(jié)果與Crank-Nicolson格式進行比較,驗證了本文中方法的有效性.

本文由3部分組成,第一部分是擴散方程的二級四階隱式Runge-Kutta方法的構(gòu)造,第二部分討論格式的穩(wěn)定性,第三部分給出具體的數(shù)值算例,并將結(jié)果與Crank-Nicolson格式和準確值進行比較,最后給出文章的結(jié)論.

1 本文中格式的構(gòu)造

1.1二級四階隱式Runge-Kutta對常微分方程組

(2)

的二級四階隱式Runge-Kutta格式為(3)式

(3)

二級四階隱式Runge-Kutta法是A-穩(wěn)定的[3].

1.2兩種差分格式的構(gòu)造方程(1)式對x變量分別用中心差分和四階緊致格式離散,對t變量保持不變,

(4)

(5)

把(4～5)式代入方程(1)式將得到如下常微分方程組(6～7)式

(6)

(7)

矩陣B嚴格對角占優(yōu)陣,因此非奇異.

常微分方程組(6～7)式用二級四階隱式Runge-Kutta格式(3)式得到本文中的格式(8～9)式

(8)

(9)

其中M=B-1C,格式(8)式為精度O(τ4+h2)格式,而(9)式為精度O(τ4+h4)格式.

下面討論格式(8)和(9)的穩(wěn)定性.

2 穩(wěn)定性分析

證明格式(8～9)式的穩(wěn)定性之前,先引入兩個引理:

引理1[2]:若A是一個N階三對角陣,

其中a,b,c是實數(shù),bc>0,則A的特征值為

由引理1直接得出矩陣A的特征值為

由此推出λA≤0.

引理2假設(shè)λi(i=1,2,…,N-1)為矩陣B-1C的特征值,xi(i=12,…,N-1)為其相應(yīng)的特征值向量,則特征值λi為實數(shù),且λi≤0.

引理2的證明下面觀察B-1C的特征值:

由此推出λi≤0.

定理1本文中差分格式(8)式和(9)式是絕對穩(wěn)定的.

定理1的證明分別記格式(8)式和(9)式的特征值為λ(8)和λ(9),則

其中λA和λi為矩陣A和B-1C的特征值.

由引理1和引理2得知λA≤0,λi≤0所以λ(8)≤1,λ(9)≤1,由此推出迭代格式(8)式和(9)式的譜半徑小于等于1,因此迭代格式(8)式和(9)式是絕對穩(wěn)定的.

3 數(shù)值實驗與結(jié)果分析

給出下面的常系數(shù)一維擴散方程初邊值問題:

該方程的準確解為u(x,t)=sin(πx)e-π2t

表1給出了當h=0.1,τ=0.1,T=1時的Crank-Nicolson格式和格式(8)式和(9)式的絕對誤差比較,最后一行給出了最大誤差,從表1可以看出格式(8)式和(9)式誤差比Crank-Nicolson格式小.

表1 h=0.1,τ=0.1,T=1時絕對誤差比較

表2給出了空間步長取不同值時的Crank-Nicolson格式和格式(8)式、(9)式的最大誤差和收斂階比較,從表中可以看出Crank-Nicolson格式和格式(8)對空間變量二階收斂,而格式(9)對空間變量四階收斂.

表2 當τ=0.001,T=1時的不同空間步長的收斂階的比較

表3給出了時間步長取不同值時Crank-Nicolson格式和格式(8)式、(9)式的最大誤差和收斂階比較,從表3可以看出格式(8)式、(9)式對時間變量四階收斂.

表3 當h=0.005,T=1時的不同空間步長的收斂階的比較

圖1給出了h=0.1,τ=0.1,T=1時Crank-Nicolson格式、格式(8)式、(9)式與準確解進行比較,從圖1中可以看出本文中格式(8)式和(9)式比Crank-Nicolson格式更接近于準確解.

4 結(jié)論

對空間變量分別采用二階中心差分格式和四階緊致差分格式進行離散話,對時間變量采用二級四階隱式Runge-Kutta公式,對擴散方程構(gòu)造了截斷誤差分別為O(h2+τ4)和O(h4+τ4)的兩種兩層絕對穩(wěn)定的隱式差分格式.試驗結(jié)果表明理論上的截斷誤差和數(shù)值試驗結(jié)果相符合.表1～3和圖1的數(shù)值試驗結(jié)果說明格式(8)式、(9)式與Crank-Nicolson格式相比,更接近于準確解.

圖1 h=0.1,τ=0.1,T=1時的3種格式與準確解比較

[1] 陸金甫,關(guān)治.偏微分方程數(shù)值解法[M].2版.北京:清華大學(xué)出版社,2004.

[2] Smith G D. Numerical solution of partial differential equations(finite difference methods)[M]. Third edition. Oxford: Oxford University Press,1996.

[3] 李瑞霞,何志慶.微分方程數(shù)值方法[M].廣州:華南理工大學(xué)出版社,2005:38-43.

[4] 開依沙爾·熱合曼.一種一維擴散方程三階精度的半離散隱式差分格式[J].甘肅聯(lián)合大學(xué)學(xué)報,2009,23(1):33-36.

[5] 開依沙爾·熱合曼.求解一維熱傳導(dǎo)方程的一種半離散差分格式[J].新疆師范大學(xué)學(xué)報,2007,26(3):142-46.

[6] 田振夫.非齊次熱傳導(dǎo)方程的高精度隱式格式[J].寧夏大學(xué)學(xué)報,1996,17(3):34-38.

[7] 葛永斌,田振夫.求解擴散方程的一種高精度隱式差分方法[J].上海理工大學(xué)學(xué)報,2005,27(2):107-111.

[8] 阿不都熱西提·阿不都外力.修正局部Crank-Nicolson法對于二維熱傳導(dǎo)方程的應(yīng)用[J].計算數(shù)學(xué),1997,19(3):267-276.

[9] Abuduwali A, Sakakihara M, Niki H. A local Crank-Nicolson method for solving the heat equation[J]. Hiroshima Mathematical Journal,1994,24(1):1-13.

[10] Sallam S, Anwar M N, Abdel-Aziz M R. Unconditionally stableC1-cubic spline collocation method for solving parabolic equations[J]. International Journal of Computer Mathematics,2004,81(7):813-821.

[11] Mohebbi A, Dehghan M. High-order compact solution of the one-dimensional heat and advection-diffusion equations[J]. Applied Mathematical Modelling,2010,34(10):3071-3084.

(責(zé)任編輯 趙燕)

Twostagefourth-orderRunge-Kuttamethodforsolvingdiffusionequation

KAYSAR·Rahman, NURMAMAT·Helil

(School of Mathematics and System Sciences, Xinjiang University, Urumqi 830046, China)

We applieed central finite difference approximation of second order and compact finite difference approximation of fourth order for discrediting spatial derivatives, and used two stage fourth order Runge-Kutta method in time direction derived two unconditionally stable implicit schemes in which local truncation error wasO(τ4+h2) andO(τ4+h4), then discussed its stability. Numerical experiment was compared with Crank-Nicolson scheme. Numerical experiment results showed that it was an efficient method for solving diffusion equation.

diffusion equation; compact scheme; two stage fourth-order Runge-Kutta method; two-level implicit scheme; Crank-Nicolson method

2014-02-13

國家自然科學(xué)基金項目(11261057)和新疆維吾爾自治區(qū)教育廳高校科研計劃重點項目(XJEDU2012I01)資助

開依沙爾·熱合曼(1978-),男,副教授,E-mail:kaysar2001@sina.com;努爾買買提·黑力力,通信作者,副教授,E-mail:nurmamat@gmail.com

1000-2375(2014)05-0476-05

O241.8

A

10.3969/j.issn.1000-2375.2014.05.019