常微分方程初值問題數(shù)值解法

會(huì)計(jì)學(xué)1常微分方程初值問題數(shù)值解法2023/1/182

在高等數(shù)學(xué)中，對(duì)于常微分方程的求解，給出了一些典型方程求解析解的基本方法，如可分離變量法、常系數(shù)齊次線性方程的解法、常系數(shù)非齊次線性方程的解法等。但能求解的常微分方程仍然是有限的，大多數(shù)的常微分方程是不能給出解析解。譬如

這個(gè)一階微分方程就不能用初等函數(shù)來表達(dá)它的解。

第1頁/共81頁2023/1/183再如，方程

的解,雖然有表可查,但對(duì)于表上沒有給出的值,仍需插值方法來計(jì)算第2頁/共81頁2023/1/184從實(shí)際問題當(dāng)中歸納出來的微分方程，通常主要依靠數(shù)值解法來解決。本章主要討論一階常微分方程初值問題

(9.1)

在區(qū)間a≤x≤b上的數(shù)值解法。

可以證明,如果函數(shù)在帶形區(qū)域R:{a≤x≤b,-∞＜y＜∞}內(nèi)連續(xù)，且關(guān)于y滿足李普希茲(Lipschitz)條件，即存在常數(shù)L(它與x,y無關(guān))使

對(duì)R內(nèi)任意x及兩個(gè)都成立,則方程(9.1)的解在a,b上存在且唯一。

第3頁/共81頁2023/1/185數(shù)值方法的基本思想

對(duì)常微分方程初值問題(9.1)式的數(shù)值解法，就是要算出精確解y(x)在區(qū)間a,b上的一系列離散節(jié)點(diǎn)處的函數(shù)值的近似值相鄰兩個(gè)節(jié)點(diǎn)的間距稱為步長，步長可以相等，也可以不等。本章總是假定h為定數(shù)，稱為定步長，這時(shí)節(jié)點(diǎn)可表示為數(shù)值解法需要把連續(xù)性的問題加以離散化，從而求出離散節(jié)點(diǎn)的數(shù)值解。

第4頁/共81頁2023/1/186

對(duì)常微分方程數(shù)值解法的基本出發(fā)點(diǎn)就是離散化。其數(shù)值解法有兩個(gè)基本特點(diǎn)，它們都采用“步進(jìn)式”，即求解過程順著節(jié)點(diǎn)排列的次序一步一步地向前推進(jìn)，描述這類算法，要求給出用已知信息計(jì)算的遞推公式。建立這類遞推公式的基本方法是在這些節(jié)點(diǎn)上用數(shù)值積分、數(shù)值微分、泰勒展開等離散化方法，對(duì)初值問題中的導(dǎo)數(shù)進(jìn)行不同的離散化處理。

第5頁/共81頁2023/1/187對(duì)于初值問題的數(shù)值解法，首先要解決的問題就是如何對(duì)微分方程進(jìn)行離散化，建立求數(shù)值解的遞推公式。遞推公式通常有兩類，一類是計(jì)算yi+1時(shí)只用到xi+1,xi和yi，即前一步的值，因此有了初值以后就可以逐步往下計(jì)算，此類方法稱為單步法；其代表是龍格—庫塔法。另一類是計(jì)算yi+1時(shí)，除用到xi+1,xi和yi以外，還要用到，即前面k步的值，此類方法稱為多步法；其代表是亞當(dāng)斯法。

第6頁/共81頁2023/1/188§9.2簡單的數(shù)值方法與基本概念9.2.1Euler公式歐拉（Euler）方法是解初值問題的最簡單的數(shù)值方法。初值問題的解y=y(x)代表通過點(diǎn)的一條稱之為微分方程的積分曲線。積分曲線上每一點(diǎn)的切線的斜率等于函數(shù)在這點(diǎn)的值。

第7頁/共81頁2023/1/189Euler法的求解過程是:從初始點(diǎn)P0(即點(diǎn)(x0,y0))出發(fā),作積分曲線y=y(x)在P0點(diǎn)上切線(其斜率為

),與x=x1直線相交于P1點(diǎn)(即點(diǎn)(x1,y1),得到y(tǒng)1作為y(x1)的近似值,如上圖所示。過點(diǎn)(x0,y0),以f(x0,y0)為斜率的切線方程為

當(dāng)時(shí),得這樣就獲得了P1點(diǎn)的坐標(biāo)。

第8頁/共81頁2023/1/1810同樣,過點(diǎn)P1(x1,y1),作積分曲線y=y(x)的切線交直線x=x2于P2點(diǎn),切線的斜率=直線方程為當(dāng)時(shí),得第9頁/共81頁2023/1/1811當(dāng)時(shí),得由此獲得了P2的坐標(biāo)。重復(fù)以上過程,就可獲得一系列的點(diǎn):P1,P1,…,Pn。對(duì)已求得點(diǎn)以=為斜率作直線

取第10頁/共81頁2023/1/1812從圖形上看,就獲得了一條近似于曲線y=y(x)的折線。這樣,從x0逐個(gè)算出對(duì)應(yīng)的數(shù)值解第11頁/共81頁2023/1/1813通常取(常數(shù)),則Euler法的計(jì)算格式

i=0,1,…,n(9.2)還可用數(shù)值微分、數(shù)值積分法和泰勒展開法推導(dǎo)Euler格式。以數(shù)值積分為例進(jìn)行推導(dǎo)。將方程的兩端在區(qū)間上積分得，

選擇不同的計(jì)算方法計(jì)算上式的積分項(xiàng)

,就會(huì)得到不同的計(jì)算公式。

(9.3)第12頁/共81頁2023/1/1814

用左矩形方法計(jì)算積分項(xiàng)

代入(9.3)式,并用yi近似代替式中y(xi)即可得到向前歐拉（Euler）公式

由于數(shù)值積分的矩形方法精度很低，所以歐拉（Euler）公式當(dāng)然很粗糙。

第13頁/共81頁2023/1/1815例9.1用歐拉法解初值問題

取步長h=0.2,計(jì)算過程保留4位小數(shù)

解:h=0.2,歐拉迭代格式

當(dāng)k=0,x1=0.2時(shí)，已知x0=0，y0=1，有

y(0.2)y1=0.2×1(4－0×1)＝0.8當(dāng)k=1,x2=0.4時(shí)，已知x1=0.2，y1=0.8，有

y(0.4)y2

=0.2×0.8×(4－0.2×0.8)＝0.6144當(dāng)k=2,x3=0.6時(shí)，已知x2=0.4，y2=0.6144，有

y(0.6)y3=0.2×0.6144×(4-0.4×0.6144)=0.4613第14頁/共81頁2023/1/1816clear;y=1,x=0,%初始化forn=1:10y=1.1*y-0.2*x/y,x=x+0.1,endy=1x=0y=1.1000x=0.1000y=1.1918x=0.2000y=1.2774x=0.3000y=1.3582x=0.4000y=1.4351x=0.5000y=1.5090x=0.6000y=1.5803x=0.7000y=1.6498x=0.8000y=1.7178x=0.9000y=1.7848x=1.0000第15頁/共81頁2023/1/18179.2.2梯形公式為了提高精度,對(duì)方程的兩端在區(qū)間上積分得，改用梯形方法計(jì)算其積分項(xiàng)，即

(9.4)

代入(7.4)式,并用近似代替式中即可得到梯形公式(9.5)

由于數(shù)值積分的梯形公式比矩形公式的精度高，因此梯形公式（9.5）比歐拉公式(9.2)的精度高一個(gè)數(shù)值方法。第16頁/共81頁2023/1/1818(9.5)

(9.5)式的右端含有未知的yi+1,它是一個(gè)關(guān)于yi+1的函數(shù)方程,這類數(shù)值方法稱為隱式方法。相反地,歐拉法是關(guān)于yi+1的一個(gè)直接的計(jì)算公式，這類數(shù)值方法稱為顯式方法。

第17頁/共81頁2023/1/18199.2.3兩步歐拉公式

對(duì)方程的兩端在區(qū)間上積分得

(9.6)

改用中矩形公式計(jì)算其積分項(xiàng)，即

代入上式,并用yi近似代替式中y(xi)即可得到兩步歐拉公式

(9.7)第18頁/共81頁2023/1/1820前面介紹過的數(shù)值方法,無論是歐拉方法,還是梯形方法，它們都是單步法,其特點(diǎn)是在計(jì)算yi+1時(shí)只用到前一步的信息yi;可是公式(7.7)中除了yi外,還用到更前一步的信息yi-1,即調(diào)用了前兩步的信息,故稱其為兩步歐拉公式

第19頁/共81頁2023/1/18219.2.4．歐拉法的局部截?cái)嗾`差衡量求解公式好壞的一個(gè)主要標(biāo)準(zhǔn)是求解公式的精度,因此引入局部截?cái)嗾`差和階數(shù)的概念。定義9.1在yi準(zhǔn)確的前提下,即時(shí),用數(shù)值方法計(jì)算yi+1的誤差,稱為該數(shù)值方法計(jì)算時(shí)yi+1的局部截?cái)嗾`差。對(duì)于歐拉公式，假定，則有而將真解y(x)在xi處按二階泰勒展開

因此有

第20頁/共81頁2023/1/1822定義9.2數(shù)值方法的局部截?cái)嗾`差為,則稱這種數(shù)值方法的階數(shù)是P。步長(h<1)越小,P越高,則局部截?cái)嗾`差越小,計(jì)算精度越高。歐拉公式的局部截?cái)嗾`差為,歐拉方法僅為一階方法。兩步歐拉公式比歐拉公式精度也是高一個(gè)數(shù)值方法,設(shè),前兩步準(zhǔn)確,則兩步歐拉公式

把y(xi-1)在xi處展開成Taylor級(jí)數(shù)，即

第21頁/共81頁2023/1/1823由

再將y(xi+1)在xi+1處進(jìn)行泰勒展開(9.8)(9.9)所以,由(9.8)和(9.9)可得兩步歐拉公式的局部截?cái)嗾`差為,即第22頁/共81頁2023/1/18249.2.5改進(jìn)的歐拉公式顯式歐拉公式計(jì)算工作量小,但精度低。梯形公式雖提高了精度,但為隱式公式,需用迭代法求解,計(jì)算工作量大。綜合歐拉公式和梯形公式便可得到改進(jìn)的歐拉公式。先用歐拉公式(9.2)求出一個(gè)初步的近似值,稱為預(yù)測值,它的精度不高,再用梯形公式(9.5)對(duì)它校正一次,即迭代一次,求得yi+1,稱為校正值,這種預(yù)測-校正方法稱為改進(jìn)的歐拉公式：（9.10）

預(yù)測

校正第23頁/共81頁2023/1/1825可以證明,公式(9.10)的精度為二階。這是一種一步顯式格式,它可以表示為嵌套形式。（9.11）或者表示成下列平均化形式（9.12）

第24頁/共81頁2023/1/18269.2.6改進(jìn)歐拉法算法實(shí)現(xiàn)（1）計(jì)算步驟①輸入,h,N②使用以下改進(jìn)歐拉法公式進(jìn)行計(jì)算

③

輸出，并使轉(zhuǎn)到

②

直至n>N結(jié)束。

第25頁/共81頁2023/1/1827（2）改進(jìn)歐拉法的流程圖

第26頁/共81頁2023/1/1828(3)

程序?qū)崿F(xiàn)(改進(jìn)歐拉法計(jì)算常微分方程初值問題

)例9.2用改進(jìn)歐拉法解初值問題區(qū)間為0,1,取步長h=0.1

解:改進(jìn)歐拉法的具體形式本題的精確解為,第27頁/共81頁2023/1/1829clearx=0,yn=1%初始化forn=1:10yp=yn+0.1*(yn-2*x/yn);%預(yù)測x=x+0.1;yc=yn+0.1*(yp-2*x/yp);yn=(yp+yc)/2%校正end第28頁/共81頁2023/1/1830例9.3對(duì)初值問題

證明用梯形公式求得的近似解為

并證明當(dāng)步長h0時(shí),yn收斂于精確解證明:解初值問題的梯形公式為∵

∴

整理成顯式

反復(fù)迭代,得到∵

∴

第29頁/共81頁2023/1/1831由于，有

∴

證畢第30頁/共81頁2023/1/1832§9.3龍格-庫塔（Runge-Kutta）法9.3.1龍格-庫塔(Runge-Kutta)法的基本思想

Euler公式可改寫成則yi+1的表達(dá)式y(tǒng)(xi+1)與的Taylor展開式的前兩項(xiàng)完全相同,即局部截?cái)嗾`差為。改進(jìn)的Euler公式又可改寫成第31頁/共81頁2023/1/1833上述兩組公式在形式上有一個(gè)共同點(diǎn):都是用f(x,y)在某些點(diǎn)上值的線性組合得出y(xi+1)的近似值yi+1,而且增加計(jì)算的次數(shù)f(x,y)的次數(shù),可提高截?cái)嗾`差的階。如歐拉公式:每步計(jì)算一次f(x,y)的值,為一階方法。改進(jìn)歐拉公式需計(jì)算兩次f(x,y)的值，它是二階方法。它的局部截?cái)嗾`差為。第32頁/共81頁2023/1/1834于是可考慮用函數(shù)f(x,y)在若干點(diǎn)上的函數(shù)值的線性組合來構(gòu)造近似公式，構(gòu)造時(shí)要求近似公式在(xi,yi)處的Taylor展開式與解y(x)在xi處的Taylor展開式的前面幾項(xiàng)重合，從而使近似公式達(dá)到所需要的階數(shù)。既避免求偏導(dǎo),又提高了計(jì)算方法精度的階數(shù)?；蛘哒f,在這一步內(nèi)多預(yù)報(bào)幾個(gè)點(diǎn)的斜率值，然后將其加權(quán)平均作為平均斜率，則可構(gòu)造出更高精度的計(jì)算格式，這就是龍格—庫塔（Runge-Kutta）法的基本思想。

第33頁/共81頁2023/1/18359.3.2二階龍格—庫塔法在上取兩點(diǎn)xi和,以該兩點(diǎn)處的斜率值k1和k2的加權(quán)平均(或稱為線性組合)來求取平均斜率k*的近似值K，即

式中:k1為xi點(diǎn)處的切線斜率值，

k2為點(diǎn)處的切線斜率值,比照改進(jìn)的歐拉法,將視為，即可得對(duì)常微分方程初值問題(9.1)式的解y=y(x),根據(jù)微分中值定理，存在點(diǎn)，使得第34頁/共81頁2023/1/1836式中K可看作是y=y(x)在區(qū)間上的平均斜率。所以可得計(jì)算公式為：（9.14）

將y(xi)在x=xi處進(jìn)行二階Taylor展開：

（9.15）

也即（9.13）第35頁/共81頁2023/1/1837將在x=xi處進(jìn)行一階Taylor展開：

將以上結(jié)果代入（9.14）得：（9.16）

對(duì)式(9.15)和(9.16)進(jìn)行比較系數(shù)后可知,只要（9.17）

成立,格式(9.14)的局部截?cái)嗾`差就等于有2階精度第36頁/共81頁2023/1/1838式(9.17)中具有三個(gè)未知量,但只有兩個(gè)方程,因而有無窮多解。若取,則p=1，這是無窮多解中的一個(gè)解，將以上所解的值代入式(9.14)并改寫可得

不難發(fā)現(xiàn)，上面的格式就是改進(jìn)的歐拉格式。凡滿足條件式（9.17）有一簇形如上式的計(jì)算格式，這些格式統(tǒng)稱為二階龍格—庫塔格式。因此改進(jìn)的歐拉格式是眾多的二階龍格—庫塔法中的一種特殊格式。第37頁/共81頁2023/1/1839若取,則，此時(shí)二階龍格-庫塔法的計(jì)算公式為

此計(jì)算公式稱為變形的二階龍格—庫塔法。式中為區(qū)間的中點(diǎn)。第38頁/共81頁2023/1/18409.3.3三階龍格-庫塔法

為了進(jìn)一步提高精度，設(shè)除外再增加一點(diǎn)并用三個(gè)點(diǎn),,的斜率k1,k2,k3加權(quán)平均得出平均斜率k*的近似值，這時(shí)計(jì)算格式具有形式:（9.18）

為了預(yù)報(bào)點(diǎn)的斜率值k3,在區(qū)間內(nèi)有兩個(gè)斜率值k1和k2可以用,可將k1,k2加權(quán)平均得出上的平均斜率,從而得到的預(yù)報(bào)值第39頁/共81頁2023/1/1841于是可得

運(yùn)用Taylor展開方法選擇參數(shù),可以使格式(9.18)的局部截?cái)嗾`差為,即具有三階精度，這類格式統(tǒng)稱為三階龍格—庫塔方法。下列是其中的一種，稱為庫塔（Kutta）公式。

（9.19）

第40頁/共81頁2023/1/1842第41頁/共81頁2023/1/18439.3.4四階龍格—庫塔法

如果需要再提高精度，用類似上述的處理方法，只需在區(qū)間上用四個(gè)點(diǎn)處的斜率加權(quán)平均作為平均斜率k*的近似值，構(gòu)成一系列四階龍格—庫塔公式。具有四階精度，即局部截?cái)嗾`差是。由于推導(dǎo)復(fù)雜，這里從略，只介紹最常用的一種四階經(jīng)典龍格—庫塔公式。

（9.20）

第42頁/共81頁2023/1/18449.3.5四階龍格—庫塔法算法實(shí)現(xiàn)(1)

計(jì)算步驟①輸入,h,N②使用龍格—庫塔公式（9.20）計(jì)算出y1③輸出，并使轉(zhuǎn)到②直至n>N結(jié)束。第43頁/共81頁2023/1/1845(2)四階龍格—庫塔算法流程圖第44頁/共81頁2023/1/1846程序?qū)崿F(xiàn)(四階龍格-庫塔法計(jì)算常微分方程初值問題)

例9.4取步長h=0.2，用經(jīng)典格式求解初值問題解:

由四階龍格-庫塔公式可得

第45頁/共81頁2023/1/1847可同樣進(jìn)行其余yi的計(jì)算。本例方程的解為，從表中看到所求的數(shù)值解具有4位有效數(shù)字。

龍格—庫塔方法的推導(dǎo)基于Taylor展開方法，因而它要求所求的解具有較好的光滑性。如果解的光滑性差，那么，使用四階龍格—庫塔方法求得的數(shù)值解，其精度可能反而不如改進(jìn)的歐拉方法。在實(shí)際計(jì)算時(shí)，應(yīng)當(dāng)針對(duì)問題的具體特點(diǎn)選擇合適的算法。

第46頁/共81頁2023/1/18489.3.6變步長的龍格-庫塔法在微分方程的數(shù)值解中，選擇適當(dāng)?shù)牟介L是非常重要的。單從每一步看，步長越小，截?cái)嗾`差就越??；但隨著步長的縮小，在一定的求解區(qū)間內(nèi)所要完成的步數(shù)就增加了。這樣會(huì)引起計(jì)算量的增大，并且會(huì)引起舍入誤差的大量積累與傳播。因此微分方程數(shù)值解法也有選擇步長的問題。以經(jīng)典的四階龍格-庫塔法(9.20)為例。從節(jié)點(diǎn)xi出發(fā)，先以h為步長求出一個(gè)近似值，記為，由于局部截?cái)嗾`差為，故有

當(dāng)h值不大時(shí)，式中的系數(shù)c可近似地看作為常數(shù)。第47頁/共81頁2023/1/1849然后將步長折半,即以為步長,從節(jié)點(diǎn)xi出發(fā),跨兩步到節(jié)點(diǎn)xi+1,再求得一個(gè)近似值,每跨一步的截?cái)嗾`差是,因此有這樣由此可得這表明以作為的近似值，其誤差可用步長折半前后兩次計(jì)算結(jié)果的偏差來判斷所選步長是否適當(dāng)?shù)?8頁/共81頁2023/1/1850當(dāng)要求的數(shù)值精度為ε時(shí)：（1）如果Δ>ε，反復(fù)將步長折半進(jìn)行計(jì)算，直至Δ<ε為止,并取其最后一次步長的計(jì)算結(jié)果作為（2）如果Δ<ε，反復(fù)將步長加倍，直到Δ>ε為止，并以上一次步長的計(jì)算結(jié)果作為。這種通過步長加倍或折半來處理步長的方法稱為變步長法。表面上看，為了選擇步長，每一步都要反復(fù)判斷Δ，增加了計(jì)算工作量，但在方程的解y(x)變化劇烈的情況下，總的計(jì)算工作量得到減少，結(jié)果還是合算的。第49頁/共81頁2023/1/1851

§9.4算法的穩(wěn)定性及收斂性9.4.1穩(wěn)定性穩(wěn)定性在微分方程的數(shù)值解法中是一個(gè)非常重要的問題。因?yàn)槲⒎址匠坛踔祮栴}的數(shù)值方法是用差分格式進(jìn)行計(jì)算的，而在差分方程的求解過程中，存在著各種計(jì)算誤差，這些計(jì)算誤差如舍入誤差等引起的擾動(dòng)，在傳播過程中，可能會(huì)大量積累，對(duì)計(jì)算結(jié)果的準(zhǔn)確性將產(chǎn)生影響。這就涉及到算法穩(wěn)定性問題。

第50頁/共81頁2023/1/1852

當(dāng)在某節(jié)點(diǎn)上xi的yi值有大小為δ的擾動(dòng)時(shí)，如果在其后的各節(jié)點(diǎn)上的值yi產(chǎn)生的偏差都不大于δ，則稱這種方法是穩(wěn)定的。

穩(wěn)定性不僅與算法有關(guān)，而且與方程中函數(shù)f(x,y)也有關(guān)，討論起來比較復(fù)雜。為簡單起見，通常只針對(duì)模型方程來討論。一般方程若局部線性化，也可化為上述形式。模型方程相對(duì)比較簡單，若一個(gè)數(shù)值方法對(duì)模型方程是穩(wěn)定的，并不能保證該方法對(duì)任何方程都穩(wěn)定，但若某方法對(duì)模型方程都不穩(wěn)定，也就很難用于其他方程的求解。第51頁/共81頁2023/1/1853先考察顯式Euler方法的穩(wěn)定性。模型方程的Euler公式為將上式反復(fù)遞推后，可得或式中第52頁/共81頁2023/1/1854

要使yi有界，其充要條件為即由于，故有（9.26）可見，如欲保證算法的穩(wěn)定，顯式Euler格式的步長h的選取要受到式（9.26）的限制。的絕對(duì)值越大，則限制的h值就越小。用隱式Euler格式，對(duì)模型方程的計(jì)算公式為，可化為第53頁/共81頁2023/1/1855由于,則恒有,故恒有因此，隱式Euler格式是絕對(duì)穩(wěn)定的（無條件穩(wěn)定）的（對(duì)任何h>0）。9.4.2收斂性常微分方程初值問題的求解，是將微分方程轉(zhuǎn)化為差分方程來求解，并用計(jì)算值yi來近似替代y(xi)，這種近似替代是否合理，還須看分割區(qū)間的長度h越來越小時(shí)，即時(shí)，是否成立。若成立，則稱該方法是收斂的，否則稱為不收斂。這里仍以Euler方法為例，來分析其收斂性。Euler格式第54頁/共81頁2023/1/1856設(shè)表示取時(shí),按Euler公式的計(jì)算結(jié)果,即Euler方法局部截?cái)嗾`差為：設(shè)有常數(shù)，則（9.27）

總體截?cái)嗾`差（9.28）

又由于f(x,y)關(guān)于y滿足李普希茨條件,即第55頁/共81頁2023/1/1857代入(9.28)上式，有再利用式（9.27），式（9.28）即上式反復(fù)遞推后，可得（9.29）

第56頁/共81頁2023/1/1858設(shè)（T為常數(shù)）

因?yàn)?/p>

所以把上式代入式（9.29），得若不計(jì)初值誤差，即，則有（9.30）

式(9.30)說明,當(dāng)時(shí),,即,所以Euler方法是收斂的，且其收斂速度為，即具有一階收斂速度。同時(shí)還說明Euler方法的整體截?cái)嗾`差為，因此算法的精度為一階。第57頁/共81頁2023/1/1859§9.5亞當(dāng)姆斯方法9.5.1亞當(dāng)姆斯格式龍格-庫塔方法是一類重要算法，但這類算法在每一步都需要先預(yù)報(bào)幾個(gè)點(diǎn)上的斜率值，計(jì)算量比較大?？紤]到計(jì)算yi+1之前已得出一系列節(jié)點(diǎn)上的斜率值，能否利用這些已知值來減少計(jì)算量呢？這就是亞當(dāng)姆斯（Adams）方法的設(shè)計(jì)思想。

第58頁/共81頁2023/1/1860設(shè)用xi,xi+1兩點(diǎn)的斜率值加權(quán)平均作為區(qū)間上的平均斜率，有計(jì)算格式

（9.21）

選取參數(shù)λ，使格式（9.21）具有二階精度。第59頁/共81頁2023/1/1861將在xi處Taylor展開

代入計(jì)算格式(9.21)化簡,并假設(shè),因此有與y(xi+1)在xi處的Taylor展開式相比較,需取

才使格式(9.21)具有二階精度。這樣導(dǎo)出的計(jì)算格式稱之為二階亞當(dāng)姆斯格式。類似地可以導(dǎo)出三階亞當(dāng)姆斯格式。

第60頁/共81頁2023/1/1862和四階亞當(dāng)姆斯格式。

(9.22)這里和下面均記上述幾種亞當(dāng)姆斯格式都是顯式的，算法比較簡單，但用節(jié)點(diǎn)的斜率值來預(yù)報(bào)區(qū)間上的平均斜率是個(gè)外推過程，效果不夠理想。為了進(jìn)一步改善精度，變外推為內(nèi)插，即增加節(jié)點(diǎn)xi+1的斜率值來得出上的平均斜率。譬如考察形如第61頁/共81頁2023/1/1863（9.23）

的隱式格式,設(shè)(9.23)右端的Taylor展開有

可見要使格式(9.23)具有二階精度,需令,這樣就可構(gòu)造二階隱式亞當(dāng)姆斯格式其實(shí)是梯形格式。類似可導(dǎo)出三階隱式亞當(dāng)姆斯格式第62頁/共81頁2023/1/1864和四階隱式亞當(dāng)姆斯格式

（9.24）

9.5.2亞當(dāng)姆斯預(yù)報(bào)-校正格式參照改進(jìn)的歐拉格式的構(gòu)造方法，以四階亞當(dāng)姆斯為例，將顯式（9.22）和隱式（9.24）相結(jié)合，用顯式公式做預(yù)報(bào)，再用隱式公式做校正，可構(gòu)成亞當(dāng)姆斯預(yù)報(bào)-校正格式（9.25）

預(yù)報(bào)：

校正：

第63頁/共81頁2023/1/1865這種預(yù)報(bào)-校正格式是四步法，它在計(jì)算yi+1時(shí)不但用到前一步的信息，而且要用到再前面三步的信息，因此它不能自行啟動(dòng)。在實(shí)際計(jì)算時(shí)，可借助于某種單步法，譬如四階龍格—庫塔法提供開始值。

第64頁/共81頁2023/1/1866例9.5取步長h=0.1,用亞當(dāng)姆斯預(yù)報(bào)-校正公式求解初值問題

的數(shù)值解。解:

用四階龍格-庫塔公式求出發(fā)值,計(jì)算得：表中的,yi和y(xi)分別為預(yù)報(bào)值、校正值和準(zhǔn)確解(),以比較計(jì)算結(jié)果的精度。再使用亞當(dāng)姆斯預(yù)報(bào)-校正公式(9.25),第65頁/共81頁2023/1/1867§9.6一階常微分方程組與高階方程我們已介紹了一階常微分方程初值問題的各種數(shù)值解法，對(duì)于一階常微分方程組，可類似得到各種解法，而高階常微分方程可轉(zhuǎn)化為一階常微分方程組來求解。

9.6.1一階常微分方程組對(duì)于一階常微分方程組的初值問題

（9.31）

可以把單個(gè)方程中的f和y看作向量來處理，這樣就可把前面介紹的各種差分算法推廣到求一階方程組初值問題中來。第66頁/共81頁2023/1/1868設(shè)為節(jié)點(diǎn)上的近似解，則有改進(jìn)的Euler格式為預(yù)報(bào)：校正：

（9.32）

又，相應(yīng)的四階龍格—庫塔格式（經(jīng)典格式）為（9.33）

第67頁/共81頁2023/1/1869式中

（9.34）

第68頁/共81頁2023/1/1870把節(jié)點(diǎn)xi上的yi和zi值代入式(7.34),依次算出

,然后把它們代入式(7.33),算出節(jié)點(diǎn)xi+1上的yi+1

和zi+1值。對(duì)于具有三個(gè)或三個(gè)以上方程的方程組的初值問題,也可用類似方法處理,只是更復(fù)雜一些而已。此外,多步方法也同樣可以應(yīng)用于求解方程組初值問題。

例7.6用改