常微分方程初值問題的數(shù)值解法中三種算法的比較

南通大學(xué)理學(xué)院 謝 舒 陳俊吉 陸海華

一、數(shù)值解法的基本思路和途徑

一階微分方程初值問題的解就是區(qū)間上的連續(xù)變量的函數(shù)，對(duì)該問題初值的解，就是求在區(qū)間[a,b]上的若干個(gè)離散點(diǎn)處的函數(shù)的近似值，比如：a≤x1≤x2＜x3＜…＜xn≤b，接著計(jì)算出解的近似值：y(x1)，y(x2)，…，y(xn)，一般取x1，x2，…，xn為等距離的點(diǎn)，即x2-x1=x3-x2=…=xn-xn-1=h或者xi=a+ih，i=1,2,3…，n，其中h稱為步長。

建立數(shù)值方法的第一步，就是把函數(shù)離散化。常用的離散化的方法有以下兩種：

（2）Taylor 展式法。在點(diǎn)x1的附近，y(x)可近似為Taylor 展開式的前P+1 項(xiàng)

其中P為一正整數(shù)。通過微分方程y'=f(x,y)，可以逐次把各階導(dǎo)數(shù)y'，y''，…，y(P)在x1處的值表示出來。

二、數(shù)值方法的導(dǎo)出與分析

1.Euler 方法

Euler 方法是最簡單的一種一階的單步法，缺點(diǎn)是精度較差，優(yōu)點(diǎn)是公式簡單，而且有比較清楚的幾何解釋，有助于直觀理解數(shù)值y1是怎樣逼近微分方程的精確解的。

特別地，當(dāng)P=1 時(shí)，近似的Taylor 多項(xiàng)式可化為上述的Euler公式。

2.梯形公式法

梯形公式（TRA）是對(duì)Euler 公式的一個(gè)改進(jìn)，梯形表達(dá)式為：

3.Runge-Kutta 方法

相對(duì)于其他常微分方程（組）的數(shù)值解法，Runge-Kutta 方法階較高，但它并沒有增加微商f(x,y)的次數(shù)，所以它相當(dāng)精確、穩(wěn)定而且容易編程。

（1）常用的二級(jí)二階Runge-Kutta 公式為：

（2）常用的三級(jí)三階Runge-Kutta 公式為

對(duì)精度要求不高的實(shí)際問題中，我們通?？梢圆捎枚?jí)二階和三級(jí)三階Runge-Kutta 方法。而對(duì)于要求較高的問題，我們通常采用如下的四級(jí)四階Runge-Kutta 方法來處理常微分方程的數(shù)值解。

（3）常用的四級(jí)四階Runge-Kutta 公式為：

上式在實(shí)際的計(jì)算過程中是應(yīng)用最廣泛的Runge-Kutta 法公式，也被稱為經(jīng)典的四階顯示Runge-Kutta 公式。

三、三種方法的實(shí)例比較

我們通過兩個(gè)具體的實(shí)例分析解析解和三種算法得到的數(shù)值解之間的關(guān)系。

1.一階常微分方程實(shí)例（1）

（1）解析解

利用y'+p(x)y=Q(x)的通解公式

可直接得到本例的通解y=Cex-x-1。將初值條件y(0)=1 代入可得本例的解為y=2ex-x-1。

（2）利用三種算法求其數(shù)值解。分別選取步長為h=0.1，利用三種算法來求解該方程在x=0.1，0.2，0.3，0.4，0.5,0.6,0.7,0.8,0.9,1.0處的數(shù)值解，將數(shù)值解的結(jié)果與解析式所求解的結(jié)果作比較。首先分別利用Euler 公式、梯形公式公式、四級(jí)四階Runge-Kutta 公式，使用Matlab 求解方程的數(shù)值解。

表1 三種解法的數(shù)值解和解析解

圖2 數(shù)值解與解析解數(shù)值比較

從圖1 和圖2 中可以很清晰地看出四級(jí)四階Runge-Kutta 方法比梯形公式法和Euler 公式法與解析解更加接近。近似度最差的是Euler 法，Euler 法的數(shù)值解與解析解的數(shù)值差最大，梯形公式法優(yōu)于Euler 法，效果最好的還是Runge-Kutta 法。Runge-Kutta 法得到的數(shù)值解明顯接近解析解的值，且隨著數(shù)值增加，梯形公式法和Euler 法與解析解的誤差越來越大，但是Runge-Kutta 法的計(jì)算過于復(fù)雜，而Euler 相對(duì)而言簡單許多。同時(shí)，我們對(duì)相關(guān)數(shù)值解和解析解進(jìn)行誤差分析，如圖3。

圖3 數(shù)值解和解析解的誤差

綜上，Runge-Kutta 法的數(shù)值解明顯精確度更高，更加接近于解析解，誤差也接近于0，但是其計(jì)算復(fù)雜程度明顯較高，計(jì)算機(jī)編程代碼比Euler 法和梯形公式法復(fù)雜許多，所以在實(shí)際應(yīng)用時(shí)要靈活采取不同的方法計(jì)算。

2.三種方法的實(shí)例比較（2）

我們考慮如下常微分方程初值問題的解析解與數(shù)值解

選取步長為h=0.1，再利用解析法和數(shù)值解法求解該方程在x=0.1，0.2，0.3，0.4，0.5,0.6,0.7,0.8,0.9,1.0 處的數(shù)值解，根據(jù)所求結(jié)果，比較數(shù)值解的結(jié)果與解析式解的結(jié)果及兩種方法的截?cái)嗾`差和誤差絕對(duì)值。

（1）解析解

（2）利用三種算法求其數(shù)值解。分別選取步長為h=0.1，再利用解析法和數(shù)值解法求解該方程在x=0.1，0.2，0.3，0.4，0.5,0.6,0.7,0.8,0.9,1.0 處的數(shù)值解，根據(jù)所求結(jié)果，比較數(shù)值解的結(jié)果與解析式解的結(jié)果及兩種方法的截?cái)嗾`差和誤差絕對(duì)值。同例1，我們可以得到該常微分方程數(shù)值解和解析解數(shù)值的比較圖以及誤差分析圖（如圖4，圖5）。

圖4 數(shù)值解與解析解數(shù)值比較

圖5 數(shù)值解和解析解的誤差

從圖4 和圖5 中我們可以看出，當(dāng)x的取值接近初值時(shí)，三個(gè)數(shù)值解與解析解非常接近，當(dāng)隨著x的逐漸增大，Euler 法的數(shù)值越來越偏離解析值，梯形法次之，Runge-Kutta法的數(shù)值與解析值最接近，并且誤差也幾乎為0。這一結(jié)果和例1 得到的結(jié)果是統(tǒng)一的。