多種微分方程數值計算方法分析

王建強，沈愉樂

(1.武漢大學測繪學院，湖北武漢 430072； 2.吳江市建設房產測量有限公司，江蘇吳江 215200)

多種微分方程數值計算方法分析

王建強1?，沈愉樂2

(1.武漢大學測繪學院，湖北武漢 430072； 2.吳江市建設房產測量有限公司，江蘇吳江 215200)

數值微分方程的數值解法是計算方法、數值分析理論中非常重要的內容，數值微分方法也是解決實際計算問題的重要方法。本文對幾種常用的數值微分方法進行了簡要的分析，并用這幾種方法對具有光滑性質的被積函數進行數值計算，龍格－庫塔方法和4階阿達姆斯方法的數值計算穩定性和計算精度都比較好。

微分方程；計算方法；數值分析；數值實驗

1 引 言

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

2 數值計算方法

本文直接給出這幾種方法的公式，具體推導過程見文獻[3]。

歐拉方法是常微分方程初值問題解中最簡單的方法。歐拉折線公式具有1階代數精度及其改進形式:

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

改進的歐拉方法，即歐拉方法的隱式公式:

改進的歐拉方法是一個預報－校正的公式，它的穩定性要有余兩步歐拉方法。龍格－庫塔方法有多種形式，并且有多種階數，常用的是標準4階龍格－庫塔法。龍格－庫塔方法的推導基于泰勒(Taylor)級數展開，它要求方程解具有良好的光滑性質。反之，如果解的光滑性差，使用4階龍格－庫塔方法求得的數值解，其精度可能不如2階改進的歐拉方法。

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

如果只計算一次，這是一個預報－校正方法。實際計算時可以進行多次迭代，迭代次數不宜過多，本文迭代計算兩次。在線性多步法的應用中，目前最常用的方法為4階阿達姆斯(Adams)外推于內插法所形成的預報－校正方法。

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

3 計算分析

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

它的解析式為:

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

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

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

第二個微分方程:

它的解析表達式為:

通過數值計算結果如圖3、圖4所示。圖3和圖4分別是步長h＝0.1和h＝0.05的計算結果，圖中的不同曲線和圖1、圖2中曲線代表的計算方法相同。計算方法的精度和前一個試驗相當，在這個試驗中，h＝0.1時，龍格－庫塔方法的計算精度比4階阿達姆斯的差，而在h＝0.05是，4階阿達姆斯的計算精度比龍格－庫塔法、辛普生方法精度都要低，計算精度的提高最不明顯。

對這兩個數值試驗，從圖1～圖4中可以看出，對于解比較光滑的曲線，龍格－庫塔法、4階阿達姆斯方法和辛普生方法的計算精度比歐拉方法、兩步歐拉方法和改進的歐拉方法要高一個數量級。辛普生方法的積分誤差曲線存在較小的波動，數值穩定性不好，4階阿達姆斯方法和龍格－庫塔法的誤差曲線隨著步長的增加而增加，數值穩定性較好。從這兩個試驗中，隨著步長的減小，計算精度的提高隨著微分方程、計算方法的不同而不同。

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

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

4 結 論

通過理論分析和兩個數值試驗可以得出，數值微分方程的隱式解比同階的顯式解的數值穩定性要好，龍格－庫塔法、辛普生方法和4階阿達姆斯方法具有較高的代數精度，可以得到比較好的結果。數值微分方程的數值解有很多種方法，建議使用計算精度較高，穩定性較好的數值計算方法，比如龍格－庫塔法、4階阿達姆斯方法等，并且有必要使用多種方法做檢核。

[1]奚梅成.數值分析方法[M].北京:中國科學技術大學出版社，2003

[2]王能超.計算方法簡明教程[M].北京:高等教育出版社，2004

[3]甄西豐.實用數值計算方法[M].北京:清華大學出版社，2006

[4]吳勃英.數值分析[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空間常微分方程初值問題弱解的一個逼近定理[J].西南大學學報(自然科學版)，30(4):49～52.關于初值問題的解的探討.

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

王建強(1981—)，男，博士研究生，研究方向為全球重力場模型研究。