第7章常微分方程數值解法7.1引言7.2簡單的數值方法與基本概念（歐拉法）7.3-庫塔方法7.5方程組和高階方程7.1引言

科學技術中常常需要求解常微分方程的定解問題.這類問題最簡單的形式，是本章將著重考察的一階方程的初值問題

我們知道，只有f(x,y)適當光滑—譬如關于y滿足萊布尼茲(Lipschitz)條件理論上就可以保證初值問題的解y＝f(x)存在并且唯一.

雖然求解常微分方程有各種各樣的解析方法，但解析方法只能用來求解一些特殊類型的方程，實際問題中歸結出來的微分方程主要靠數值解法.

所謂數值解法,就是尋求解y(x)在一系列離散節點上的近似值y1,y2,,yn,yn+1,.相鄰兩個節點的間距hn=xn+1-xn稱為步長.今后如不特別說明，總是假定hi=h(i=1,2,)為常數,這時節點為xn=x0+nh(i=0,1,2,)(等距節點).截去最后一項,可得這就是著名的(顯式)歐拉(Euler)公式.若初值y0已知，則依公式(2.1)可逐次逐步算出各點數值解.即7.2簡單的數值方法與基本概念7.2.1歐拉法與后退歐拉法用向前差商代替導數

例1

用歐拉公式求解初值問題

解取步長h=0.1，歐拉公式的具體形式為其中xn=nh=0.1n(n=0,1,,10),已知y0=1,由此式可得依次計算下去，部分計算結果見下表.與準確解相比，可看出歐拉公式的計算結果精度很差.

xn

歐拉公式數值解yn準確解y(xn)

誤差0.20.40.60.81.01.1918181.3582131.5089661.6497831.7847701.1832161.3416411.4832401.6124521.7320510.0086020.0165720.0257260.0373310.052719

歐拉公式具有明顯的幾何意義,就是用折線近似代替方程的解曲線，因而常稱公式(2.1)為歐拉折線法.

還可以通過幾何直觀來考察歐拉方法的精度.假設yn=y(xn),即頂點Pn落在積分曲線y=y(x)上，那么，按歐拉方法做出的折線PnPn+1便是y=y(x)過點Pn的切線.從圖形上看,這樣定出的頂點Pn+1顯著地偏離了原來的積分曲線，可見歐拉方法是相當粗糙的.

為了分析計算公式的精度，通常可用泰勒展開將y(xn+1)在xn處展開，則有在yn=y(xn)的前提下，f(xn,yn)=f(xn,y(xn))=y(xn).于是可得歐拉法(2.1)的公式誤差為稱為此方法的局部截斷誤差.截去最后一項,可得這就是著名的(隱式)后退歐拉(Euler)公式.若初值y0已知，則依公式(2.5)可逐次逐步算出各點數值解.即用向后差商代替導數

后退的歐拉公式與歐拉公式有著本質的區別,后者是關于yn+1的一個直接計算公式，這類公式稱作是顯式的；前者公式的右端含有未知的yn+1，它實際上是關于yn+1的一個函數方程,這類方程稱作是隱式的.

顯式與隱式兩類方法各有特點，考了到數值穩定性等其他因素，人們有時需要選用隱式方法，但使用顯式算法遠比隱式方便.

隱式方程通常用迭代法求解，而迭代過程的實質是逐步顯式化.7.2.2梯形方法

為得到比歐拉法精度高的計算公式，將歐拉公式和后退的歐拉公式進行算數平均，可得稱為梯形公式.

梯形公式是隱式單步法，用迭代法求解，同后退的歐拉方法一樣，仍用歐拉法提供迭代初值。7.2.3單步法的局部截斷誤差與階

初值問題(1.1),(1.2)的單步法可用一般形式表示為其中多元函數與f(x,y

)有關，當含有yn+1時，方法是隱式的，若不含yn+1則為顯式方法，所以顯式單步法可表示為(x,y,h)稱為增量函數，例如對歐拉法(2.1)有它的局部截斷誤差(2.3)已由給出,對一般顯式單步法則可如下定義.

定義1

設y(x)是初值問題(1.1),(1.2)的準確解,稱為顯式單步法(2.10)的局部截斷誤差.

Tn+1之所以稱為局部的，是假設在xn前各步沒有誤差.當yn=y(xn)時，計算一步，則有

所以，局部截斷誤差可理解為用方法(2.10)計算一步的誤差，也即公式(2.10)中用準確解y(x)代替數值解產生的公式誤差.根據定義,顯然歐拉法的局部截斷誤差為顯然Tn+1=O(h2).一般情形的定義如下

定義2

設y(x)是初值問題的準確解，若存在最大整數p使顯式單步法(2.10)的局部截斷誤差滿足則稱方法(2.10)具有p階精度.若將(2.10)展開式寫成則稱為局部截斷誤差主項.

以上定義對隱式單步法(2.9)也是適用的.例如，對后退歐拉法(2.5)其局部截斷誤差為這里p=1是1階方法，局部截斷誤差主項為同樣對梯形法(2.7)有所以梯形方法(2.7)是二階的.其局部截斷誤差主項為7.2.4改進的歐拉公式我們看到，梯形方法雖然提高了精度，但其算法復雜，在應用迭代公式(2.9)進行實際計算時，每迭代一次，都要重新計算函數f(x,y

)的值，而迭代又要反復進行若干次，計算量很大，而且往往難以預測.為了控制計算量，通常只迭代一兩次就轉入下一步的計算，這就簡化了算法.具體地說，我們先用歐拉公式求得一個初步的近似值，稱之為預測值，此預測值的精度可能很差，再用梯形公式(2.7)將它校正一次，即按(2.8)式迭代一次，這個結果稱之為校正值.這樣建立的預測—校正系統通常稱為改進的歐拉公式:或表示為下列平均化形式(2.13)預測校正

例2

用改進的歐拉法解例1中的初值問題(2.2).

解仍取步長h=0.1,改進的歐拉公式為部分計算結果見下表

xnyn

誤差

xnyn

誤差00.20.41.1.1840961.34336000.0000880.0017190.60.81.01.4859561.6164761.7378690.0027160.0040240.005818同例1中的歐拉法的計算結果比較，明顯改善了精度.

xn

歐拉公式

改進歐拉公式

yn

誤差

yn

誤差0.00.20.40.60.81.01.1.1918181.3582131.5089661.6497831.784770

00.0086020.016572

0.0257260.0373310.05271911.1840961.3433601.4859561.6164761.737869

00.0000880.0017190.0027160.0040240.0058187.3龍格—庫塔方法

對許多實際問題來說，歐拉公式與改進歐拉公式精度還不能滿足要求，為此從另一個角度來分析這兩個公式的特點，從而探索一條構造高精度方法的途徑.

7.3.1顯式龍格—庫塔法的一般形式

上節給出了顯式單步法的表達式(2.10),其局部截斷誤差為O(hp+1)，對歐拉法Tn+1=O(h2)，即方法為p=1階，若用改進歐拉法(2.13)，它可表為此時增量函數為一般說來，點數r越多，精度越高，上式右端相當于增量函數(xn,yn,h)，為得到便于計算的顯式方法，可類似于改進歐拉法(3.1),(3.2)，將公式表示為其中這里ci,ai,bij均為常數.(3.4)和(3.5)稱為r級顯式龍格-庫塔(Runge-Kutta)法,簡稱R-K方法.7.3.2二階顯式R-K方法

對r=2的R-K方法，由(3.4),(3.5)式可得如下計算公式這里c1,c2,a2,b21均為待定常數，我們希望適當選取這些系數，使公式階數p盡量高.

則由此可以看出在改進的歐拉公式中相當于取(xn,yn),(xn+1,yn+1)兩點處斜率的平均值，近似代替平均斜率，其精度比歐拉公式提高了.c1=c2=1/2,a2=b21=a=1時，即為改進的歐拉公式：稱為變形的歐拉公式,也可以表示為:

如取a=1，則c1=0,c2=1,a2=b21=1/2.得計算公式

對r=2的R-K公式能否使局部誤差提高到O(h4)?為此需把k2多展開一項，展開式中的項是不能通過選擇參數削掉的，實際上要使h3

的項為零，需增加3個方程，要確定4個參數c1,c2,λ2,μ21，這是不可能的.故r＝2的顯式R-K方法的階只能是p=2，而不能得到三階公式.

四階R-K方法的每一步需要計算四次函數值f，可以證明其局部截斷誤差為O(h5).

7.3.3標準四階顯式R-K方法

然而值得指出的是，龍格-庫塔方法的推導基于泰勒展開方法，因而它要求所求的解具有較好的光滑性質.反之，如果解的光滑性差，那么，使用龍格-庫塔方法求得的數值解，其精度可能反而不如改進的歐拉方法.實際計算時，我們應當針對問題的具體特點選擇合適的算法.7.3.4變步長的龍格-庫塔方法

單從每一步看，步長越小，截斷誤差就越小，但隨著步長的縮小，在一定求解范圍內所要完成的步數就增加了.步數的增加不但引起計算量的增大，而且可能導致舍入誤差的嚴重積累.因此同積分的數值計算一樣，微分方程的數值解法也有個選擇步長的問題.

在選擇步長時，需要考慮兩個問題：

1.怎樣衡量和檢驗計算結果的精度？

2.如何依據所獲得的精度處理步長？我們考察四階R-K公式(3.13)，從節點xn出發，先以h為步長求出一個近似值，記為，由于公式的局部截斷誤差為O(h5)，故然后將步長折半，即取為步長,從xn跨兩步到xn+1，再求得一個近似值，每跨一步的局部截斷誤差是

，因此有比較(3.14)式和(3.15)式我們看到，步長折半后，誤差大約減少到1/16，即有由此易得下列事后估計式這樣，我們可以通過檢查步長，折半前后兩次計算結果的偏差來判定所選的步長是否合適，具體地說，將區分以下兩種情況處理：這種通過加倍或折半處理步長的方法稱為變步長方法.表面上看，為了選擇步長，每一步的計算量增加了，但總體考慮往往是合算的.

1.對于給定的精度ε,如果Δ>ε,我們反復將步長折半計算,直至Δ<ε為止,這時取最終得到的作為結果；

2.如果Δ<ε，我們將反復將步長作加倍計算，直至Δ>ε為止，這時再將步長折半計算一次，就得到所要的結果.7.5方程組和高階方程7.5.1一階方程組

前面我們研究了單個方程y=f

的數值解，只要把y和f理解為向量，那么，所提供的各種計算公式即可應用到一階方程組的情形.

考察一階方程組的初值問題，初始條件給為若采用向量的記號，記(向量)則上述方程組的初值問題可表示為求解這一初值問題的四階龍格-庫塔公式為(向量)式中(向量)或表示為(分量)其中(分量)這里yin是第i個因變量yi(x)在節點xn=x0+nh的近似值.為了幫助理解這一公式的計算過程，我們再考察兩個方程的特殊情形這時四階龍格-庫塔公式具有形式其中