王建強(qiáng)，沈愉樂(lè)

(1.武漢大學(xué)測(cè)繪學(xué)院，湖北武漢 430072； 2.吳江市建設(shè)房產(chǎn)測(cè)量有限公司，江蘇吳江 215200)

多種微分方程數(shù)值計(jì)算方法分析

王建強(qiáng)1?，沈愉樂(lè)2

(1.武漢大學(xué)測(cè)繪學(xué)院，湖北武漢 430072； 2.吳江市建設(shè)房產(chǎn)測(cè)量有限公司，江蘇吳江 215200)

數(shù)值微分方程的數(shù)值解法是計(jì)算方法、數(shù)值分析理論中非常重要的內(nèi)容，數(shù)值微分方法也是解決實(shí)際計(jì)算問(wèn)題的重要方法。本文對(duì)幾種常用的數(shù)值微分方法進(jìn)行了簡(jiǎn)要的分析，并用這幾種方法對(duì)具有光滑性質(zhì)的被積函數(shù)進(jìn)行數(shù)值計(jì)算，龍格－庫(kù)塔方法和4階阿達(dá)姆斯方法的數(shù)值計(jì)算穩(wěn)定性和計(jì)算精度都比較好。

微分方程；計(jì)算方法；數(shù)值分析；數(shù)值實(shí)驗(yàn)

1 引 言

在科學(xué)研究和工程應(yīng)用中，所建立的數(shù)學(xué)模型多是常微分方程或微分方程組，但是除了少數(shù)特殊類型的微分方程能用解析方法求得其精確解外，大多數(shù)情況下要得出解的解析表達(dá)式是極其困難的，因此，就需要用數(shù)值逼近方法求得其近似解。微分方程組的數(shù)值解的存在性和穩(wěn)定性取決于被積函數(shù)的特性和初值。求初值問(wèn)題的數(shù)值解法可區(qū)分為兩大類:單步法和多步法。常用的方法中歐拉(Euler)方法、改進(jìn)的歐拉方法、龍格－庫(kù)塔方法(Runge－Kutta)等[1～4]是單步法的典型代表，線性多步法是多步法的典型代表，對(duì)于一些特別的數(shù)值微分方程使用這些方法效果很差[5，6]。微分方程的數(shù)值解法有顯式解和隱式解法，一般來(lái)說(shuō)，隱式解要優(yōu)于顯式解[4]。歐拉方法是一種最簡(jiǎn)單的單步法，計(jì)算量小，但精度比較低。一般的初值問(wèn)題，多采用改進(jìn)的歐拉方法，因?yàn)樗臄?shù)值穩(wěn)定性和計(jì)算精度都比一般的歐拉方法好。龍格－庫(kù)塔方法是一類應(yīng)用較廣的高精度單步法，當(dāng)解充分光滑時(shí)的4階龍格－庫(kù)塔方法一般可以達(dá)到很高的精度。常微分方程的初值對(duì)計(jì)算方法的收斂是有影響的[7]。為了更好地比較這幾種常用的方法，本文采用這幾種數(shù)值方法對(duì)被積函數(shù)光滑連續(xù)，初值精確的微分方程做了數(shù)值試驗(yàn)。

2 數(shù)值計(jì)算方法

本文直接給出這幾種方法的公式，具體推導(dǎo)過(guò)程見(jiàn)文獻(xiàn)[3]。

歐拉方法是常微分方程初值問(wèn)題解中最簡(jiǎn)單的方法。歐拉折線公式具有1階代數(shù)精度及其改進(jìn)形式:

式(1)中，k為節(jié)點(diǎn)序列，h為步長(zhǎng)，f(x，y)為已知函數(shù)。兩步歐拉方法具有2階代數(shù)精度:

改進(jìn)的歐拉方法，即歐拉方法的隱式公式:

改進(jìn)的歐拉方法是一個(gè)預(yù)報(bào)－校正的公式，它的穩(wěn)定性要有余兩步歐拉方法。龍格－庫(kù)塔方法有多種形式，并且有多種階數(shù)，常用的是標(biāo)準(zhǔn)4階龍格－庫(kù)塔法。龍格－庫(kù)塔方法的推導(dǎo)基于泰勒(Taylor)級(jí)數(shù)展開(kāi)，它要求方程解具有良好的光滑性質(zhì)。反之，如果解的光滑性差，使用4階龍格－庫(kù)塔方法求得的數(shù)值解，其精度可能不如2階改進(jìn)的歐拉方法。

以上4種方法都是單步法，由于線性多步法需要利用前面的若干個(gè)點(diǎn)，所以初始計(jì)算時(shí)使用單步法計(jì)算得到若干個(gè)點(diǎn)，再利用多步法進(jìn)行計(jì)算，一般使用龍格－庫(kù)塔作為初始計(jì)算方法。利用辛普森公式建立的遞推公式為:

如果只計(jì)算一次，這是一個(gè)預(yù)報(bào)－校正方法。實(shí)際計(jì)算時(shí)可以進(jìn)行多次迭代，迭代次數(shù)不宜過(guò)多，本文迭代計(jì)算兩次。在線性多步法的應(yīng)用中，目前最常用的方法為4階阿達(dá)姆斯(Adams)外推于內(nèi)插法所形成的預(yù)報(bào)－校正方法。

由于多步法計(jì)算需要前面幾步的數(shù)值解，公式(5)需要前面的3步，公式(6)需要前面的4步，因此數(shù)值計(jì)算時(shí)用龍格－庫(kù)塔作為初始計(jì)算起始時(shí)的幾步數(shù)值解。

3 計(jì)算分析

利用以上6種方法，計(jì)算了兩個(gè)微分方程，這兩個(gè)微分方程都有解析表達(dá)式，因此可以求出積分誤差。第一個(gè)微分方程為:

它的解析式為:

圖1 試驗(yàn)1誤差曲線(h＝0.1)

圖2 試驗(yàn)1誤差曲線(h＝0.05)

通過(guò)數(shù)值計(jì)算結(jié)果如圖1、圖2所示。圖1和圖2分別是步長(zhǎng)h＝0.1和h＝0.05的計(jì)算結(jié)果。在圖1和圖2中，(a)圖中的點(diǎn)線曲線是由歐拉折線公式計(jì)算的結(jié)果，精度最低，可以精確到小數(shù)點(diǎn)后1位。實(shí)線是由兩步歐拉方法計(jì)算的結(jié)果，它比歐拉折線公式得到結(jié)果精度要高，可以精確到小數(shù)點(diǎn)后2位，從圖中可以看出它的誤差有明顯的振蕩，而不是隨著k值的增加而增加，表明計(jì)算方法的數(shù)值穩(wěn)定性不好。虛線表示改進(jìn)的歐拉方法，它的計(jì)算精度和兩步歐拉方法基本相同，在有些點(diǎn)上甚至不如兩步歐拉方法，但是它的誤差隨著 k的增加而增加，數(shù)值計(jì)算的穩(wěn)定性較好。(b)圖中的點(diǎn)線是龍格－庫(kù)塔計(jì)算的結(jié)果，它的誤差隨著k的增加而增加，數(shù)值計(jì)算的穩(wěn)定性較好，實(shí)線是辛普森公式建立的遞推公式的結(jié)果，它的誤差有較小的振蕩，虛線的是4階阿達(dá)姆斯方法求得的數(shù)值解。從圖1和圖2中可以看出，當(dāng)步長(zhǎng)減小后，4階阿達(dá)姆斯方法計(jì)算精度提高得最明顯。

第二個(gè)微分方程:

它的解析表達(dá)式為:

通過(guò)數(shù)值計(jì)算結(jié)果如圖3、圖4所示。圖3和圖4分別是步長(zhǎng)h＝0.1和h＝0.05的計(jì)算結(jié)果，圖中的不同曲線和圖1、圖2中曲線代表的計(jì)算方法相同。計(jì)算方法的精度和前一個(gè)試驗(yàn)相當(dāng)，在這個(gè)試驗(yàn)中，h＝0.1時(shí)，龍格－庫(kù)塔方法的計(jì)算精度比4階阿達(dá)姆斯的差，而在h＝0.05是，4階阿達(dá)姆斯的計(jì)算精度比龍格－庫(kù)塔法、辛普生方法精度都要低，計(jì)算精度的提高最不明顯。

對(duì)這兩個(gè)數(shù)值試驗(yàn)，從圖1～圖4中可以看出，對(duì)于解比較光滑的曲線，龍格－庫(kù)塔法、4階阿達(dá)姆斯方法和辛普生方法的計(jì)算精度比歐拉方法、兩步歐拉方法和改進(jìn)的歐拉方法要高一個(gè)數(shù)量級(jí)。辛普生方法的積分誤差曲線存在較小的波動(dòng)，數(shù)值穩(wěn)定性不好，4階阿達(dá)姆斯方法和龍格－庫(kù)塔法的誤差曲線隨著步長(zhǎng)的增加而增加，數(shù)值穩(wěn)定性較好。從這兩個(gè)試驗(yàn)中，隨著步長(zhǎng)的減小，計(jì)算精度的提高隨著微分方程、計(jì)算方法的不同而不同。

圖3 試驗(yàn)2誤差曲線(h＝0.1)

圖4 試驗(yàn)2誤差曲線(h＝0.05)

4 結(jié) 論

通過(guò)理論分析和兩個(gè)數(shù)值試驗(yàn)可以得出，數(shù)值微分方程的隱式解比同階的顯式解的數(shù)值穩(wěn)定性要好，龍格－庫(kù)塔法、辛普生方法和4階阿達(dá)姆斯方法具有較高的代數(shù)精度，可以得到比較好的結(jié)果。數(shù)值微分方程的數(shù)值解有很多種方法，建議使用計(jì)算精度較高，穩(wěn)定性較好的數(shù)值計(jì)算方法，比如龍格－庫(kù)塔法、4階阿達(dá)姆斯方法等，并且有必要使用多種方法做檢核。

[1]奚梅成.數(shù)值分析方法[M].北京:中國(guó)科學(xué)技術(shù)大學(xué)出版社，2003

[2]王能超.計(jì)算方法簡(jiǎn)明教程[M].北京:高等教育出版社，2004

[3]甄西豐.實(shí)用數(shù)值計(jì)算方法[M].北京:清華大學(xué)出版社，2006

[4]吳勃英.數(shù)值分析[M].北京:高等教育出版社，2007

[5]Iserles A.On the Method of Neumann Series for Highly Oscillating Equations[J].BIT，2004，44:473～488

[6]Iserles A.On the Global Error of Discretization Methods for Highly－Oscillatory Ordinary Differential Equations[J]. BIT，2002，42:561～599

[7]陳清明.Banach空間常微分方程初值問(wèn)題弱解的一個(gè)逼近定理[J].西南大學(xué)學(xué)報(bào)(自然科學(xué)版)，30(4):49～52.關(guān)于初值問(wèn)題的解的探討.

Numerical Calculation Analysis on Different Methods of Differential Equation

Wang JianQiang1，Shen YuLe2
(1.School of Geodesy and Geomatics，Wuhan University，Wuhan 430079，China；2.Wujiang construction and real estate of survey Co.，Ltd.Wujiang 215200，China)

Numerical solution of differential equations is the very important content of numerical calculation method and numerical analysis，and numerical differentiation method is important to solve practical computing problems.In this paper，several commonly used methods of numerical differentiation are given a brief analysis and use these methods with the integrand is smooth to carry out numerical experiments，Runge－Kutta method and fourth－order Adams method Numerical calculation are better than others with the stability and calculation accuracy.

differential equation；calculation method；numerical analysis；numerical experiment

1672－8262(2010)04－117－03

O175

A

2009—12—26

王建強(qiáng)(1981—)，男，博士研究生，研究方向?yàn)槿蛑亓?chǎng)模型研究。