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

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

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

我們知道微分方程是關(guān)于未知函數(shù)及導(dǎo)數(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)是一個關(guān)于

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

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

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

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

點所求的計算解，稱

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

為該數(shù)值方法在

點的局部截斷誤差。

顯然局部截斷誤差是從

到

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

準確的前提下，

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

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

計算

、再從

計算

、…、直到從

計算

這

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

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

之前的

是準確時，局部截斷誤差

與總體截斷誤差

是相同的。

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

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

定義7.3

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

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

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

定義7.4

設(shè)

是初值問題在

點的兩個近似初值，

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

，當(dāng)步長滿足

且

時，有下式成立

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

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

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

開始求解的近似解序列

會落入積分曲線

一個鄰近通道內(nèi)。如下圖所示：初值問題

定理7.1設(shè)

在區(qū)域

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

對任意

和任意

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

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

1.Euler公式

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

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

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

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

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

表示，如圖。

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

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

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

初值條件

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

以曲線在

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

得到

，即

歐拉方法

這里

稱為步長。同樣以

點處的

值為方向，前進到直線

到

，即

如圖所示，再從

得到

，從

得到

，…

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

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

處用向前差商推導(dǎo)而得； 2)Euler公式的本質(zhì)是用當(dāng)前點的斜率來預(yù)測曲線下一個點

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

例7.1

用Euler公式求解初值問題

在

處的近似值，取步長

。

解：

由Euler公式

詳細計算結(jié)果列表如下

該問題的解析解為

，因此

在

處，

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

2.截斷誤差

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

設(shè)

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

則

為Euler公式在

處的局部截斷誤差，考慮

在

處的Taylor展開式

將

帶入此式并注意到

是

的解析解，可得歐拉方法

從而局部截斷誤差為

其中

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

同價。2)總體截斷誤差

設(shè)

是由Euler公式得到的，稱

為Euler公式在

處的總體截斷誤差，記作

若

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

有

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

證明：設(shè)

是Euler公式在

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

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

考慮第二部分，由Lipschitiz條件

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

注意到

，當(dāng)

時有

再注意到

則

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

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

3.后退Euler公式

前述討論的Euler公式是對

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

由此導(dǎo)出另一種形式的Euler公式，稱為后退Euler公式

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

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

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

，它實際上是關(guān)于

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

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

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

一個點的信息來構(gòu)造

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

兩個點的信息來構(gòu)造

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

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

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

設(shè)

是初值問題的解，因此有

，由此得到

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

稱

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

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

設(shè)

是初值問題的準確解，記

為梯形公式在

處的局部截斷誤差，考慮

在

處的Taylor展開式

從第二個Taylor展開式可得到

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

其局部截斷誤差為

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

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

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

3)改進Euler公式

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

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

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

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

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

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

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

解：

預(yù)測

校正

詳細計算結(jié)果列表如下

n10.021.020000001.0204000020.041.041208001.0416160830.061.063248401.0636647240.081.086138021.0865627550.101.109894011.11032732歐拉方法

該問題的解析解為

，因此

在

處，

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

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

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

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

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

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

令

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

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

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

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

1.Runge—Kutta法的基本思想

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

和

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

考慮在

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

令

，上式可寫成如下形式

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

如果把

看作是

相對于

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

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

作為校正量

的近似值。

改進Euler公式可以表示為

而改進Euler方法是用兩個點

、

處的信息構(gòu)造

和

的算術(shù)平均值作為

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

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

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

點處的值

作為預(yù)測值代入

計算出

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

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

內(nèi)多預(yù)測幾個點的

值，并用它們的加權(quán)平均來作為

的近似值，即

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

2．二階Runge—Kutta法

考慮用

和

的加權(quán)平均值來近似

的一般形式，設(shè)

要使其局部截斷誤差達到

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

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

取

，就是改進Euler公式。

2)Heun公式

取

，則得到Heun公式。

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

取

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

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

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

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

對應(yīng)C點，真正的積分曲線在

點的準確值

對應(yīng)A點，中點公式的計算值

對應(yīng)B點附近。

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

其中

就是圖中三角形與

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

更接近于準確的

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

對應(yīng)的B點更接近于準確解的A點。

對于積分曲線在

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

3.四階Runge—Kutta法

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

的加權(quán)平均值來近似

，因此令

欲使局部截斷誤差達到

，即

，類似于二階Runge—Kutta法，經(jīng)過一系列較為復(fù)雜的推導(dǎo)過程，得到有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ù)值解也有一個選擇步長問題?？墒褂萌缦路椒ǎ?/p>

設(shè)

是從

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

是從

點以

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

，如果相鄰兩次計算滿足

則停止計算，以最后計算結(jié)果為

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

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

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

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

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

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

利用前面得到的多步信息

來計算

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

其中

為常數(shù)，

為

的近似值，

，當(dāng)

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

的線性組合來構(gòu)造

。(7.26)

定義7.2

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

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

設(shè)有如下3個線性二步法公式

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

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

的信息來構(gòu)造計算

的公式，其形式如下

我們首先給出

的(7.28)式(當(dāng)

時是Euler公式)，注意到

，將

在

點Taylor展開

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

如果我們將

在

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

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

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

則局部截斷誤差

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

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

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

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