求解擴(kuò)散方程的二級(jí)四階隱式Runge-Kutta方法

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

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

求解擴(kuò)散方程的二級(jí)四階隱式Runge-Kutta方法

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

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

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

擴(kuò)散方程;緊致格式;二級(jí)四階Runge-Kutta方法;兩層隱格式;Crank-Nicolson格式

0 引言

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

(1)

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

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

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

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

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

1.1二級(jí)四階隱式Runge-Kutta對(duì)常微分方程組

(2)

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

(3)

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

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

(4)

(5)

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

(6)

(7)

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

常微分方程組(6～7)式用二級(jí)四階隱式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)定性之前,先引入兩個(gè)引理:

引理1[2]:若A是一個(gè)N階三對(duì)角陣,

其中a,b,c是實(shí)數(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í)數(shù),且λi≤0.

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

由此推出λi≤0.

定理1本文中差分格式(8)式和(9)式是絕對(duì)穩(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)式是絕對(duì)穩(wěn)定的.

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

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

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

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

表1 h=0.1,τ=0.1,T=1時(shí)絕對(duì)誤差比較

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

表2 當(dāng)τ=0.001,T=1時(shí)的不同空間步長的收斂階的比較

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

表3 當(dāng)h=0.005,T=1時(shí)的不同空間步長的收斂階的比較

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

4 結(jié)論

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

圖1 h=0.1,τ=0.1,T=1時(shí)的3種格式與準(zhǔn)確解比較

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

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

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

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

[8] 阿不都熱西提·阿不都外力.修正局部Crank-Nicolson法對(duì)于二維熱傳導(dǎo)方程的應(yīng)用[J].計(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é)基金項(xiàng)目(11261057)和新疆維吾爾自治區(qū)教育廳高?？蒲杏?jì)劃重點(diǎn)項(xiàng)目(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