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

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

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

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

常微分方程的初值問(wèn)題

考慮如下常微分方程初值問(wèn)題

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

在點(diǎn)

處的近似值

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

，稱為步長(zhǎng)，并令

，從

求得

，再?gòu)?/p>

求得

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

各點(diǎn)的數(shù)值解

。

(7.1)(7.2)

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

的方法僅使用

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

的方法需要

則稱為該方法求解初值問(wèn)題的多步法。

一般的單步法形如

在(7.3-1)式中右端無(wú)

項(xiàng)可直接計(jì)算，稱為顯式單步法；而(7.3-2)公式無(wú)法直接計(jì)算

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

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

例如

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

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

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

初值問(wèn)題(7.4)(7.4-1)(7.4-2)

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

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

點(diǎn)所求的計(jì)算解，稱

為該數(shù)值方法的整體截?cái)嗾`差或總體截?cái)嗾`差；稱

為該數(shù)值方法在

點(diǎn)的局部截?cái)嗾`差。

顯然局部截?cái)嗾`差是從

到

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

準(zhǔn)確的前提下，

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

所謂總體截?cái)嗾`差是算法從

計(jì)算

、再?gòu)?/p>

計(jì)算

、…、直到從

計(jì)算

這

n+1個(gè)步驟所積累的總誤差。

因此局部截?cái)嗾`差與總體截?cái)嗾`差是不同的但有關(guān)聯(lián)，特別如果在計(jì)算

之前的

是準(zhǔn)確時(shí)，局部截?cái)嗾`差

與總體截?cái)嗾`差

是相同的。

求解微分方程初值問(wèn)題數(shù)值方法的收斂性問(wèn)題是假定每一步計(jì)算都沒(méi)有舍入誤差的前提下，討論步長(zhǎng)

時(shí)，方法的總體截?cái)嗾`差是否趨于零的問(wèn)題。初值問(wèn)題

定義7.3

如果數(shù)值方法的局部截?cái)嗾`差可表示為

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

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

定義7.4

設(shè)

是初值問(wèn)題在

點(diǎn)的兩個(gè)近似初值，

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

，當(dāng)步長(zhǎng)滿足

且

時(shí)，有下式成立

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

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

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

開始求解的近似解序列

會(huì)落入積分曲線

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

定理7.1設(shè)

在區(qū)域

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

對(duì)任意

和任意

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

此處列出定理7.1的原因是因?yàn)橐院蟮暮芏嘟Y(jié)論均是以Lipschitz條件為前提的。初值問(wèn)題(7.7)7.2Euler方法

1.Euler公式

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

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

所表示的積分曲線上每一點(diǎn)

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

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

表示，如圖。

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

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

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

初值條件

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

以曲線在

處的切線為方向，前進(jìn)到直線

得到

，即

歐拉方法

這里

稱為步長(zhǎng)。同樣以

點(diǎn)處的

值為方向，前進(jìn)到直線

到

，即

如圖所示，再?gòu)?/p>

得到

，從

得到

，…

按此過(guò)程我們得到如下迭代公式：

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

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

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

例7.1

用Euler公式求解初值問(wèn)題

在

處的近似值，取步長(zhǎng)

。

解：

由Euler公式

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

該問(wèn)題的解析解為

，因此

在

處，

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

2.截?cái)嗾`差

現(xiàn)在我們對(duì)Euler方法進(jìn)行誤差分析，也就是討論Euler公式的截?cái)嗾`差。1)局部截?cái)嗾`差

設(shè)

是微分方程初值問(wèn)題的準(zhǔn)確解，令

則

為Euler公式在

處的局部截?cái)嗾`差，考慮

在

處的Taylor展開式

將

帶入此式并注意到

是

的解析解，可得歐拉方法

從而局部截?cái)嗾`差為

其中

稱為誤差的主項(xiàng)，從(7.9)式可以看出局部截?cái)嗾`差與

同價(jià)。2)總體截?cái)嗾`差

設(shè)

是由Euler公式得到的，稱

為Euler公式在

處的總體截?cái)嗾`差，記作

若

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

有

則Euler公式的總體截?cái)嗾`差為歐拉方法(7.9)(7.10)

證明：設(shè)

是Euler公式在

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

其中第一部分是局部截?cái)嗾`差有

考慮第二部分，由Lipschitiz條件

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

注意到

，當(dāng)

時(shí)有

再注意到

則

有界，即總體截?cái)嗾`差與h同階，也就說(shuō)明

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

3.后退Euler公式

前述討論的Euler公式是對(duì)

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

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

后退Euler公式與(向前)Euler公式的局部截?cái)嗾`差是同階的，都是

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

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

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

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

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

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

一個(gè)點(diǎn)的信息來(lái)構(gòu)造

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

兩個(gè)點(diǎn)的信息來(lái)構(gòu)造

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

4.梯形公式和改進(jìn)Euler公式1)梯形公式

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

設(shè)

是初值問(wèn)題的解，因此有

，由此得到

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

稱

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

2)梯形公式的局部截?cái)嗾`差

設(shè)

是初值問(wèn)題的準(zhǔn)確解，記

為梯形公式在

處的局部截?cái)嗾`差，考慮

在

處的Taylor展開式

從第二個(gè)Taylor展開式可得到

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

其局部截?cái)嗾`差為

從而局部截?cái)嗾`差與h3同階，其誤差主項(xiàng)為

仿照Euler方法總體截?cái)嗾`差的推導(dǎo)過(guò)程，可以證明梯形法的總體截?cái)嗾`差為

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

3)改進(jìn)Euler公式

梯形公式的局部截?cái)嗾`差要比Euler公式高一階，提高了精度，但因?yàn)樘菪喂绞请[式的，必須通過(guò)求解方程得方式得到

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

先用Euler公式求得一個(gè)初步的近似值

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

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

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

例7.2用改進(jìn)Euler公式求解例7.1中的初值問(wèn)題。

解：

預(yù)測(cè)

校正

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

n10.021.020000001.0204000020.041.041208001.0416160830.061.063248401.0636647240.081.086138021.0865627550.101.109894011.11032732歐拉方法

該問(wèn)題的解析解為

，因此

在

處，

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

4)改進(jìn)Euler公式的局部截?cái)嗾`差

雖然改進(jìn)Euler公式是簡(jiǎn)化的梯形公式，僅用一次迭代作為所求隱式方程的解，但是所得算法的局部截?cái)嗾`差卻與

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

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

考慮改進(jìn)Euler公式一步產(chǎn)生的誤差，即局部截?cái)嗾`差，

令

是梯形公式的局部截?cái)嗾`差，

是改進(jìn)Euler公式的局部截?cái)嗾`差，則歐拉方法

即改進(jìn)Euler公式局部截?cái)嗾`差與

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

1.Runge—Kutta法的基本思想

前節(jié)給出的Euler公式和改進(jìn)Euler公式的局部截?cái)嗾`差分別與

和

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

考慮在

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

令

，上式可寫成如下形式

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

如果把

看作是

相對(duì)于

的校正量(通常是未知的)，則Euler方法實(shí)際上是用一個(gè)點(diǎn)

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

作為校正量

的近似值。

改進(jìn)Euler公式可以表示為

而改進(jìn)Euler方法是用兩個(gè)點(diǎn)

、

處的信息構(gòu)造

和

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

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

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

再仔細(xì)觀察，在改進(jìn)Euler公式中，把

點(diǎn)處的值

作為預(yù)測(cè)值代入

計(jì)算出

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

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

內(nèi)多預(yù)測(cè)幾個(gè)點(diǎn)的

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

的近似值，即

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

2．二階Runge—Kutta法

考慮用

和

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

的一般形式，設(shè)

要使其局部截?cái)嗾`差達(dá)到

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

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

取

，就是改進(jìn)Euler公式。

2)Heun公式

取

，則得到Heun公式。

3)修正的Euler公式或中點(diǎn)公式

取

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

我們可以通過(guò)中點(diǎn)公式的幾何意義來(lái)體會(huì)二階Euler公式為何更為有效。首先考察積分曲線在

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

從圖上看出Euler公式的計(jì)算值

對(duì)應(yīng)C點(diǎn)，真正的積分曲線在

點(diǎn)的準(zhǔn)確值

對(duì)應(yīng)A點(diǎn)，中點(diǎn)公式的計(jì)算值

對(duì)應(yīng)B點(diǎn)附近。

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

其中

就是圖中三角形與

重合的一條直角邊的邊長(zhǎng)，該邊長(zhǎng)較之于

更接近于準(zhǔn)確的

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

對(duì)應(yīng)的B點(diǎn)更接近于準(zhǔn)確解的A點(diǎn)。

對(duì)于積分曲線在

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

3.四階Runge—Kutta法

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

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

，因此令

欲使局部截?cái)嗾`差達(dá)到

，即

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

例7.3用標(biāo)準(zhǔn)四階Runge—Kutta公式求解例7.1中的初值問(wèn)題

解：

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

4．變步長(zhǎng)Runge—Kutta法

在實(shí)際計(jì)算中，步長(zhǎng)h的選擇是一個(gè)十分重要的問(wèn)題。若步長(zhǎng)h取得太小，將增加許多不必要的計(jì)算量；若步長(zhǎng)

取得太大,其計(jì)算精度很難得到保證。因此，仿照數(shù)值積分方法，常微分方程的數(shù)值解也有一個(gè)選擇步長(zhǎng)問(wèn)題?？墒褂萌缦路椒ǎ?/p>

設(shè)

是從

點(diǎn)以h為步長(zhǎng)使用四階Runge—Kutta公式計(jì)算得到的，

是從

點(diǎn)以

為步長(zhǎng)兩次使用四階Runge—Kutta法計(jì)算過(guò)程得到的。對(duì)于給定精度

，如果相鄰兩次計(jì)算滿足

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

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

所謂多步法就是在逐步迭代的求解過(guò)程中，計(jì)算

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

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

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

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

利用前面得到的多步信息

來(lái)計(jì)算

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

其中

為常數(shù)，

為

的近似值，

，當(dāng)

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

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

。(7.26)

定義7.2

設(shè)y(x)是微分方程初值問(wèn)題的準(zhǔn)確解，如果多步法公式的局部截?cái)嗾`差滿足

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

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

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

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

的信息來(lái)構(gòu)造計(jì)算

的公式，其形式如下

我們首先給出

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

時(shí)是Euler公式)，注意到

，將

在

點(diǎn)Taylor展開

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

如果我們將

在

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

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

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

則局部截?cái)嗾`差

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

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

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

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