計算方法第7微分方程數(shù)值解

計算方法北京工業(yè)大學應用數(shù)理學院楊中華第七章常微分方程數(shù)值解

有很多實際問題是用微分方程建立其數(shù)學模型，但是目前只能對少量典型的微分方程求出相應的解析形式的解，而對于絕大多數(shù)的微分方程問題，很難或者根本不可能得到它的解析解。因此，研究求解微分方程的數(shù)值解是非常有意義的。

我們知道微分方程是關于未知函數(shù)及導數(shù)的方程，包括常微分方程和偏微分方程。常微分方程所涉及的未知函數(shù)是一元函數(shù)；如果未知函數(shù)是多元的，則稱為偏微分方程，偏微分方程數(shù)值解有另外獨立的課程介紹。本章僅討論常微分方程初值問題的數(shù)值解方法。7.1

常微分方程的初值問題

考慮如下常微分方程初值問題

所謂求初值問題的數(shù)值解就是求其解析解

在點

處的近似值

一般采用等距節(jié)點的方式，也就是取

，稱為步長，并令

，從

求得

，再從

求得

，依次求解就可得到y(tǒng)=y(x)在

各點的數(shù)值解

。

(7.1)(7.2)

定義7.1如果確定數(shù)值解

的方法僅使用

則稱該方法為求解初值問題的單步法；如果確定

的方法需要

則稱為該方法求解初值問題的多步法。

一般的單步法形如

在(7.3-1)式中右端無

項可直接計算，稱為顯式單步法；而(7.3-2)公式無法直接計算

必須通過解方程得到，稱為隱式單步法。初值問題(7.3)(7.3-1)(7.3-2)如

求解微分方程的線性多步法形如

例如

這兩個計算公式均是線性二步法，而(7.4-1)式可以直接計算

，是顯式線性二步法；注意到(7.4-2)式中右端因此(7.4-2)是一個關于

的方程，是隱式線性多步法。

初值問題(7.4)(7.4-1)(7.4-2)

定義7.2設y(x)是初值問題的解析解，

是按某數(shù)值方法在

點所求的計算解，稱

為該數(shù)值方法的整體截斷誤差或總體截斷誤差；稱

為該數(shù)值方法在

點的局部截斷誤差。

顯然局部截斷誤差是從

到

一個算法步驟所產(chǎn)生的誤差，也就是在

準確的前提下，

的誤差；初值問題(7.5)(7.6)

所謂總體截斷誤差是算法從

計算

、再從

計算

、…、直到從

計算

這

n+1個步驟所積累的總誤差。

因此局部截斷誤差與總體截斷誤差是不同的但有關聯(lián)，特別如果在計算

之前的

是準確時，局部截斷誤差

與總體截斷誤差

是相同的。

求解微分方程初值問題數(shù)值方法的收斂性問題是假定每一步計算都沒有舍入誤差的前提下，討論步長

時，方法的總體截斷誤差是否趨于零的問題。初值問題

定義7.3

如果數(shù)值方法的局部截斷誤差可表示為

則稱這種數(shù)值方法的階數(shù)是p。

數(shù)值方法的階數(shù)決定了該方法的收斂性和收斂速度。

定義7.4

設

是初值問題在

點的兩個近似初值，

是某數(shù)值方法從這兩個初值開始所求的近似解，如果存在常數(shù)

，當步長滿足

且

時，有下式成立

則稱該數(shù)值方法是穩(wěn)定的。初值問題

定義7.4其含義是如果求解微分方程的數(shù)值方法以

的任意兩個初值作為起點，所得兩近似解的誤差將被兩個初值的差所界定，進而可以推出，穩(wěn)定的數(shù)值方法在求解過程中具有自動校正的性質，也就是從

開始求解的近似解序列

會落入積分曲線

一個鄰近通道內。如下圖所示：初值問題

定理7.1設

在區(qū)域

內連續(xù)且關于y滿足Lipschitz條件：存在常數(shù)L，使

對任意

和任意

成立，則初值問題(7.1)存在連續(xù)可微且唯一的解。

此處列出定理7.1的原因是因為以后的很多結論均是以Lipschitz條件為前提的。初值問題(7.7)7.2Euler方法

1.Euler公式

先從研究一階微分方程初值問題的幾何直觀入手，問題(7.1)式可以表示為平面問題：

這表明該微分方程的解析解

所表示的積分曲線上每一點

切線斜率等于函數(shù)

的值，也就是說，積分曲線上每一點的切線方向可用二元函數(shù)

表示，如圖。

求解此微分方程初值問題的困難在于二元函數(shù)

中的y是所求的未知函數(shù)，除初始點的

外其他點函數(shù)值均未知，對于其他點的函數(shù)近似值，可使用如下方法求解：

初值條件

作為積分曲線的出發(fā)點

以曲線在

處的切線為方向，前進到直線

得到

，即

歐拉方法

這里

稱為步長。同樣以

點處的

值為方向，前進到直線

到

，即

如圖所示，再從

得到

，從

得到

，…

按此過程我們得到如下迭代公式：

這就是著名的Euler公式。 1)Euler公式是在

處用向前差商推導而得； 2)Euler公式的本質是用當前點的斜率來預測曲線下一個點

的位置，然后用折線來近似所求的積分曲線； 3)當h很小時折線對積分曲線的偏離不會太大，而且Euler 公式具有自動校正的特點，即折線會向積分曲線靠攏歐拉方法(7.8)

例7.1

用Euler公式求解初值問題

在

處的近似值，取步長

。

解：

由Euler公式

詳細計算結果列表如下

該問題的解析解為

，因此

在

處，

，3位有效數(shù)字n10.021.0200000020.041.0408000030.061.0624160040.081.0848643250.101.10816161歐拉方法

2.截斷誤差

現(xiàn)在我們對Euler方法進行誤差分析，也就是討論Euler公式的截斷誤差。1)局部截斷誤差

設

是微分方程初值問題的準確解，令

則

為Euler公式在

處的局部截斷誤差，考慮

在

處的Taylor展開式

將

帶入此式并注意到

是

的解析解，可得歐拉方法

從而局部截斷誤差為

其中

稱為誤差的主項，從(7.9)式可以看出局部截斷誤差與

同價。2)總體截斷誤差

設

是由Euler公式得到的，稱

為Euler公式在

處的總體截斷誤差，記作

若

關于y滿足Lipschitiz條件，即存在正數(shù)L，使得對一切

有

則Euler公式的總體截斷誤差為歐拉方法(7.9)(7.10)

證明：設

是Euler公式在

處由準確的y(x)計算得到的結果，則有

其中第一部分是局部截斷誤差有

考慮第二部分，由Lipschitiz條件

將兩部分綜合考慮則有歐拉方法

注意到

，當

時有

再注意到

則

有界，即總體截斷誤差與h同階，也就說明

也就說明Euler方法是收斂的，且為一階的方法。歐拉方法

3.后退Euler公式

前述討論的Euler公式是對

用向前差商近似導數(shù)得到的，如果使用向后差商，則有

由此導出另一種形式的Euler公式，稱為后退Euler公式

后退Euler公式與(向前)Euler公式的局部截斷誤差是同階的，都是

，但后退Euler公式與Euler公式有著本質的區(qū)別。Euler公式是關于

的一個直接計算公式，這類公式稱為顯式；而后退Euler公式的右端含有未知的

，它實際上是關于

的一個方程，這類公式為隱式。歐拉方法(7.11)

通常選用第二章的迭代法求解該隱式方程，而只需迭代一、二步就可以達到精度要求。

向前Euler公式僅僅根據(jù)

一個點的信息來構造

，而后退Euler公式是根據(jù)

兩個點的信息來構造

，可以預期后退Euler公式將優(yōu)于向前Euler公式，實際計算表明，后退Euler公式更穩(wěn)定。

4.梯形公式和改進Euler公式1)梯形公式

前述的Euler公式是用直觀方式導出的，其實Euler公式也可以用積分方式導出，加以改進可以得到另一個更有效的求解微分方程的公式，以下是導出過程。歐拉方法

設

是初值問題的解，因此有

，由此得到

在此公式推導過程中，是用左矩形的面積近似定積分的計算，顯然定積分的計算如果使用梯形積分公式會更精確些，由此得到

稱

為求解常微分方程的梯形公式。歐拉方法(7.11)

2)梯形公式的局部截斷誤差

設

是初值問題的準確解，記

為梯形公式在

處的局部截斷誤差，考慮

在

處的Taylor展開式

從第二個Taylor展開式可得到

將其代入第一個Taylor展開式可得到歐拉方法

其局部截斷誤差為

從而局部截斷誤差與h3同階，其誤差主項為

仿照Euler方法總體截斷誤差的推導過程，可以證明梯形法的總體截斷誤差為

，即梯形法是一個二階方法。歐拉方法(7.13)

3)改進Euler公式

梯形公式的局部截斷誤差要比Euler公式高一階，提高了精度，但因為梯形公式是隱式的，必須通過求解方程得方式得到

，算法相對復雜計算量也大。為了簡化算法，通常在使用求解非線性方程的迭代法時只迭代一兩次就轉入下一步計算，簡化算法如下：

先用Euler公式求得一個初步的近似值

，稱為預測值，預測值的精度可能較差，再用梯形公式將其校正一次得

稱為校正值，迭代格式如下：

該迭代格式稱為預測—校正公式，也稱為改進Euler公式。歐拉方法(7.14)

例7.2用改進Euler公式求解例7.1中的初值問題。

解：

預測

校正

詳細計算結果列表如下

n10.021.020000001.0204000020.041.041208001.0416160830.061.063248401.0636647240.081.086138021.0865627550.101.109894011.11032732歐拉方法

該問題的解析解為

，因此

在

處，

，有5位有效數(shù)字，比Euler公式計算結果多2位有效數(shù)字。歐拉方法

4)改進Euler公式的局部截斷誤差

雖然改進Euler公式是簡化的梯形公式，僅用一次迭代作為所求隱式方程的解，但是所得算法的局部截斷誤差卻與

同階，因此改進Euler法也是二階的方法。證明如下：

設y=y(x)是初值問題的準確解，若記

考慮改進Euler公式一步產(chǎn)生的誤差，即局部截斷誤差，

令

是梯形公式的局部截斷誤差，

是改進Euler公式的局部截斷誤差，則歐拉方法

即改進Euler公式局部截斷誤差與

同階，從而可以知道改進Euler法是二階的方法。歐拉方法7.3Runge-Kutta法

1.Runge—Kutta法的基本思想

前節(jié)給出的Euler公式和改進Euler公式的局部截斷誤差分別與

和

同階，顯然兩公式所得計算解的精度后者將優(yōu)于前者，本節(jié)先探討兩者優(yōu)劣的原因，再引導我們構造出更為有效的計算公式。

考慮在

上對y=y(x)應用微分中值定理

令

，上式可寫成如下形式

則Euler公式可表示為(7.15)(7.16)

如果把

看作是

相對于

的校正量(通常是未知的)，則Euler方法實際上是用一個點

處的信息(包括步長和積分曲線在此點上的斜率)構造

作為校正量

的近似值。

改進Euler公式可以表示為

而改進Euler方法是用兩個點

、

處的信息構造

和

的算術平均值作為

的近似值，也就是說改進Euler方法用了更多的信息構造校正量

的近似值，當然如此改進的Euler方法效果應該優(yōu)于Euler方法。龍格|庫塔方法(7.17)

再仔細觀察，在改進Euler公式中，把

點處的值

作為預測值代入

計算出

，然后再用兩值的平均值去近似

。進一步，如果能在區(qū)間

內多預測幾個點的

值，并用它們的加權平均來作為

的近似值，即

可以預期會構造出具有更高精度的計算公式，這就是Runge—Kutta法的基本思想。龍格|庫塔方法(7.18)

2．二階Runge—Kutta法

考慮用

和

的加權平均值來近似

的一般形式，設

要使其局部截斷誤差達到

，則其參數(shù)必須滿足以下方程組

注意到方程組有三個方程，4個未知量，因此方程組有無窮多個解，不同的解將對應不同的計算公式，通常有如下三種參數(shù)的取法。龍格|庫塔方法(7.19) 1)改進Euler法

取

，就是改進Euler公式。

2)Heun公式

取

，則得到Heun公式。

3)修正的Euler公式或中點公式

取

，所得公式稱為修正的Euler公式或中點公式，即龍格|庫塔方法(7.20)(7.21)

我們可以通過中點公式的幾何意義來體會二階Euler公式為何更為有效。首先考察積分曲線在

點附近是凸函數(shù)的情況，如圖：龍格|庫塔方法

從圖上看出Euler公式的計算值

對應C點，真正的積分曲線在

點的準確值

對應A點，中點公式的計算值

對應B點附近。

事實上，當積分曲線是凸函數(shù)且h足夠小時，將有龍格|庫塔方法

其中

就是圖中三角形與

重合的一條直角邊的邊長，該邊長較之于

更接近于準確的

，從而經(jīng)此校正之后的計算解

對應的B點更接近于準確解的A點。

對于積分曲線在

點附近是凹函數(shù)的情況，也容易看出中點公式優(yōu)于Euler公式，如下圖所示。

3.四階Runge—Kutta法

按照Runge—Kutta法的基本思想，構造四階Runge—Kutta法也就是用

的加權平均值來近似

，因此令

欲使局部截斷誤差達到

，即

，類似于二階Runge—Kutta法，經(jīng)過一系列較為復雜的推導過程，得到有13個參數(shù)的11個方程組(見文獻[7])。由于方程的個數(shù)少于未知量的個數(shù)，因此方程組有無窮多個解，此處僅給出兩個常用的四階公式。龍格|庫塔方法(7.22) (1)標準四階Runge—Kutta公式龍格|庫塔方法(7.23) (2)Gill公式(7.24)

例7.3用標準四階Runge—Kutta公式求解例7.1中的初值問題

解：

與準確解比較高達7位有效數(shù)字。龍格|庫塔方法

4．變步長Runge—Kutta法

在實際計算中，步長h的選擇是一個十分重要的問題。若步長h取得太小，將增加許多不必要的計算量；若步長

取得太大,其計算精度很難得到保證。因此，仿照數(shù)值積分方法，常微分方程的數(shù)值解也有一個選擇步長問題。可使用如下方法：

設

是從

點以h為步長使用四階Runge—Kutta公式計算得到的，

是從

點以

為步長兩次使用四階Runge—Kutta法計算過程得到的。對于給定精度

，如果相鄰兩次計算滿足

則停止計算，以最后計算結果為

的近似值；否則步長折半繼續(xù)計算直到滿足精度要求為止。龍格|庫塔方法(7.25)7.4線性多步法

所謂多步法就是在逐步迭代的求解過程中，計算

時已經(jīng)算出近似值

，如果能夠充分利用前面多步的信息來預測

，則可能獲得較高的精度，這就是所謂的線性多步法的基本思想。

1.線性多步法的一般公式

利用前面得到的多步信息

來計算

最常用的方法是線性多步法，其一般公式為

其中

為常數(shù)，

為

的近似值，

，當

時該式為顯式，否則為隱式。多步法的本質就是用

的線性組合來構造

。(7.26)

定義7.2

設y(x)是微分方程初值問題的準確解，如果多步法公式的局部截斷誤差滿足

則稱該多步法具有p階精度，或稱方法是p階的。顯然p愈大結果愈好。 2.線性二步法

設有如下3個線性二步法公式

前兩個公式是具有2階精度的顯式線性二步法，第三個公式是利用Simpson求積公式求得具有4階精度的隱式線性二步法。線性多步法(7.27)

3.Adams外推公式 Adams外推公式利用

的信息來構造計算

的公式，其形式如下

我們首先給出

的(7.28)式(當

時是Euler公式)，注意到

，將

在

點Taylor展開

將其代入到(7.28)式，則有

如果我們將

在

點Taylor展開線性多步法(7.28)(7.29)(7.30)

欲使(7.28)式的精度足夠高，應使

中的p盡可能的大，也就是使(7.29)和(7.30)兩式相同的項數(shù)盡可能的多，顯然應取

則局部截斷誤差

，由此得到具有二階精度的Adams公式：線性多步法(7.31)

Adams外推公式是一種顯式格式,所以也稱為Adams顯式公式，仿照上述方法我們還可以推導出

情況的Adams公式，此處不再賃述(最常用的是

情況4階精度的Adams公式)。線性多步法(7.