第七章常微分方程數(shù)值解

第七章常微分方程數(shù)值解第1頁，共93頁，2023年，2月20日，星期三第七章

常微分方程的數(shù)值解法

問題提出倒葫蘆形狀容器壁上的刻度問題.對于如圖所示圓柱形狀容器壁上的容積刻度,可以利用圓柱體體積公式其中直徑D為常數(shù).由于體積V與相對于容器底部的任意高度H的函數(shù)關(guān)系明確,因此在容器上可以方便地標出容器刻度,而對于幾何形狀不是規(guī)則的容器,比如倒葫蘆形狀容器壁上如何標出刻度呢?下表是經(jīng)過測量得到部分容器高度與直徑的關(guān)系.第2頁，共93頁，2023年，2月20日，星期三x1o根據(jù)上表的數(shù)據(jù),可以擬合出倒葫蘆形狀容器的圖,建立如圖所示的坐標軸后,問題即為如何根據(jù)任意高度x標出容器體積V的刻度,由微元思想分析可知H00.20.40.60.81.0D0.040.110.260.561.041.17其中x表示高度,直徑D是高度x的函數(shù),記為D(x),因此得到如下微分方程初值問題第3頁，共93頁，2023年，2月20日，星期三第七章

常微分方程的數(shù)值解法

只要求解上述方程,就可求出體積V與高度x之間的函數(shù)關(guān)系,從而可標出容器壁上容積的刻度,但問題是(7.0)中的函數(shù)D(x)無解析表達式,我們根本無法求出其解析解.怎樣求解不能用解析方法求解的微分方程初制值問題,就是本章要討論的重點.(7.0)第4頁，共93頁，2023年，2月20日，星期三第七章

常微分方程的數(shù)值解法

7.1引言包含自變量、未知函數(shù)及未知函數(shù)的導(dǎo)數(shù)或微分的方程稱為微分方程。在微分方程中,自變量的個數(shù)只有一個,稱為常微分方程。自變量的個數(shù)為兩個或兩個以上的微分方程叫偏微分方程。微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程的階數(shù)。如果未知函數(shù)y及其各階導(dǎo)數(shù)都是一次的,則稱它是線性的,否則稱為非線性的。第5頁，共93頁，2023年，2月20日，星期三

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

這個一階微分方程就不能用初等函數(shù)及其積分來表達它的解。

第6頁，共93頁，2023年，2月20日，星期三再如，方程

的解,雖然有表可查,但對于表上沒有給出的值,仍需插值方法來計算第7頁，共93頁，2023年，2月20日，星期三從實際問題當中歸納出來的微分方程，通常主要依靠數(shù)值解法來解決。本章主要討論一階常微分方程初值問題

(7.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))使

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

第8頁，共93頁，2023年，2月20日，星期三常微分方程表示方法第9頁，共93頁，2023年，2月20日，星期三在微分方程中,自變量的個數(shù)只有一個,稱為常微分方程.自變量的個數(shù)為兩個或兩個以上的微分方程叫偏微分方程。第10頁，共93頁，2023年，2月20日，星期三例如第11頁，共93頁，2023年，2月20日，星期三當x=0時,y=1,可得c=1特解當x=0時,y=1,可得c=-1特解兩邊積分通解常微分方程解法回顧-第12頁，共93頁，2023年，2月20日，星期三例：設(shè)降落傘從跳傘塔下落后，所受空氣阻力與速度成正比，并設(shè)降落傘離開跳傘塔時（t=0）速度為零，求降落傘下落速度與時間的函數(shù)關(guān)系。解：設(shè)降落傘下落速度為v(t),降落傘在空中下落時，同時受到重力P與阻力R的作用，重力大小為mg，方向與v一致；阻力大小為kv（k為比例系數(shù)），方向與v相反，從而降落傘所受外力為F=mg-kv根據(jù)牛頓第二定律F=ma(a為加速度），得函數(shù)v(t)應(yīng)滿足的方程為

按題意，初始條件為分離變量后得兩端積分考慮到mg-kv>0即或?qū)⒊跏紬l件代如入上式得于是所求特解為第13頁，共93頁，2023年，2月20日，星期三第14頁，共93頁，2023年，2月20日，星期三7.2數(shù)值方法的基本思想對常微分方程初值問題(7.1)式的數(shù)值解法，就是要算出精確解y(x)在區(qū)間a,b上的一系列離散節(jié)點處的函數(shù)值的近似值。相鄰兩個節(jié)點的間距稱為步長，步長可以相等，也可以不等。本章總是假定h為定數(shù)，稱為定步長，這時節(jié)點可表示為數(shù)值解法需要把連續(xù)性的問題加以離散化，從而求出離散節(jié)點的數(shù)值解。

第15頁，共93頁，2023年，2月20日，星期三

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

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

第17頁，共93頁，2023年，2月20日，星期三7.3歐拉（Euler）法7.3.1Euler公式歐拉（Euler）方法是解初值問題的最簡單的數(shù)值方法。初值問題的解y=y(x)代表通過點的一條稱之為微分方程的積分曲線。積分曲線上每一點的切線的斜率等于函數(shù)在這點的值。

第18頁，共93頁，2023年，2月20日，星期三Euler法的求解過程是:從初始點P0(即點(x0,y0))出發(fā),作積分曲線y=y(x)在P0點上切線(其斜率為),與x=x1直線相交于P1點(即點(x1,y1),得到y(tǒng)1作為y(x1)的近似值,如上圖所示。過點(x0,y0),以f(x0,y0)為斜率的切線方程為

當時,得這樣就獲得了P1點的坐標。

第19頁，共93頁，2023年，2月20日，星期三同樣,過點P1(x1,y1),作積分曲線y=y(x)的切線交直線x=x2于P2點,切線的斜率=直線方程為當時,得第20頁，共93頁，2023年，2月20日，星期三當時,得由此獲得了P2的坐標。重復(fù)以上過程,就可獲得一系列的點:P1,P1,…,Pn。對已求得點以=為斜率作直線

取第21頁，共93頁，2023年，2月20日，星期三從圖形上看,就獲得了一條近似于曲線y=y(x)的折線。這樣,從x0逐個算出對應(yīng)的數(shù)值解第22頁，共93頁，2023年，2月20日，星期三通常取(常數(shù)),則Euler法的計算格式

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

選擇不同的計算方法計算上式的積分項,就會得到不同的計算公式。

(7.3)第23頁，共93頁，2023年，2月20日，星期三

公式推導(dǎo)第24頁，共93頁，2023年，2月20日，星期三

第25頁，共93頁，2023年，2月20日，星期三

用左矩形方法計算積分項

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

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

第26頁，共93頁，2023年，2月20日，星期三

第27頁，共93頁，2023年，2月20日，星期三例7.1用歐拉法解初值問題

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

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

當k=0,x1=0.2時，已知x0=0，y0=1，有y(0.2)y1=0.2×1(4－0×1)＝0.8當k=1,x2=0.4時，已知x1=0.2，y1=0.8，有y(0.4)y2

=0.2×0.8×(4－0.2×0.8)＝0.6144當k=2,x3=0.6時，已知x2=0.4，y2=0.6144，有y(0.6)y3=0.2×0.6144×(4-0.4×0.6144)=0.4613第28頁，共93頁，2023年，2月20日，星期三例7.2用Euler法求初值問題y'=x–y2

y(0)=0的數(shù)值解（取h=0.1n=5)解：∵f(x,y)=x-y2;x0=y0=0;h=0.1∴由Euler法的遞推公式得：yn+1=yn+0.1(xn–y2n)

yn=0n=0,1,2,3,4,5由上式計算所得數(shù)據(jù)n012345…xn00.10.20.30.40.5…yn000.010.029990.05999

0.0995…

第29頁，共93頁，2023年，2月20日，星期三二、隱式Euler(尤拉)格式第30頁，共93頁，2023年，2月20日，星期三7.3.2梯形公式為了提高精度,對方程的兩端在區(qū)間上積分得，改用梯形方法計算其積分項，即

(7.4)

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

由于數(shù)值積分的梯形公式比矩形公式的精度高，因此梯形公式（7.5）比歐拉公式(7.2)的精度高一個數(shù)值方法。第31頁，共93頁，2023年，2月20日，星期三(7.5)

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

第32頁，共93頁，2023年，2月20日，星期三例7.3用梯形公式求下面初值問題的解在x=0.01上的值y(0.01)y'=y

y(0)=1解：取h=0.01,由梯形公式得y=exe0.01=1.010050167第33頁，共93頁，2023年，2月20日，星期三7.3.3兩步歐拉公式

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

(7.6)

改用中矩形公式計算其積分項，即

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

(7.7)第34頁，共93頁，2023年，2月20日，星期三前面介紹過的數(shù)值方法,無論是歐拉方法,還是梯形方法，它們都是單步法,其特點是在計算yi+1時只用到前一步的信息yi;可是公式(7.7)中除了yi外,還用到更前一步的信息yi-1,即調(diào)用了前兩步的信息,故稱其為兩步歐拉公式

第35頁，共93頁，2023年，2月20日，星期三7.3.4．歐拉法的局部截斷誤差衡量求解公式好壞的一個主要標準是求解公式的精度,因此引入局部截斷誤差和階數(shù)的概念。定義7.1在yi準確的前提下,即時,用數(shù)值方法計算yi+1的誤差,稱為該數(shù)值方法計算時yi+1的局部截斷誤差。對于歐拉公式，假定，則有而將真解y(x)在xi處按二階泰勒展開

因此有

第36頁，共93頁，2023年，2月20日，星期三定義7.2數(shù)值方法的局部截斷誤差為,則稱這種數(shù)值方法的階數(shù)是P。步長(h<1)越小,P越高,則局部截斷誤差越小,計算精度越高。歐拉公式的局部截斷誤差為,歐拉方法僅為一階方法。兩步歐拉公式比歐拉公式精度也是高一個數(shù)值方法,設(shè),前兩步準確,則兩步歐拉公式

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

將在xn處按二階泰勒展開略去余項用近似值yn代替y(xn)有第37頁，共93頁，2023年，2月20日，星期三由

再將y(xi+1)在xi處進行泰勒展開(7.8)(7.9)所以,由(7.8)和(7.9)可得兩步歐拉公式的局部截斷誤差為,即第38頁，共93頁，2023年，2月20日，星期三7.3.5改進的歐拉公式顯式歐拉公式計算工作量小,但精度低。梯形公式雖提高了精度,但為隱式公式,需用迭代法求解,計算工作量大。綜合歐拉公式和梯形公式便可得到改進的歐拉公式。先用歐拉公式(7.2)求出一個初步的近似值,稱為預(yù)測值,它的精度不高,再用梯形公式(7.5)對它校正一次,即迭代一次,求得yi+1,稱為校正值,這種預(yù)測-校正方法稱為改進的歐拉公式：（7.10）

預(yù)測

校正第39頁，共93頁，2023年，2月20日，星期三可以證明,公式(7.10)的精度為二階。這是一種一步顯式格式,它可以表示為嵌套形式。（7.11）或者表示成下列平均化形式（7.12）

第40頁，共93頁，2023年，2月20日，星期三7.3.6改進歐拉法算法實現(xiàn)（1）計算步驟①輸入,h,N②使用以下改進歐拉法公式進行計算

③

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

②

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

第41頁，共93頁，2023年，2月20日，星期三（2）改進歐拉法的流程圖

第42頁，共93頁，2023年，2月20日，星期三(3)

程序?qū)崿F(xiàn)(見附錄AA-15改進歐拉法計算常微分方程初值問題)例7.2用改進歐拉法解初值問題區(qū)間為0,1,取步長h=0.1

解:改進歐拉法的具體形式本題的精確解為,計算見P158列表所示

第43頁，共93頁，2023年，2月20日，星期三例7.3對初值問題

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

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

∴

整理成顯式

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

∴

第44頁，共93頁，2023年，2月20日，星期三由于，有

∴

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

第48頁，共93頁，2023年，2月20日，星期三7.4.2二階龍格—庫塔法在上取兩點xi和,以該兩點處的斜率值k1和k2的加權(quán)平均(或稱為線性組合)來求取平均斜率k*的近似值K，即

式中:k1為xi點處的切線斜率值，k2為點處的切線斜率值,比照改進的歐拉法,將視為，即可得對常微分方程初值問題(7.1)式的解y=y(x),根據(jù)微分中值定理，存在點，使得第49頁，共93頁，2023年，2月20日，星期三式中K可看作是y=y(x)在區(qū)間上的平均斜率。所以可得計算公式為：（7.14）

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

（7.15）

也即（7.13）第50頁，共93頁，2023年，2月20日，星期三將在x=xi處進行一階Taylor展開：

將以上結(jié)果代入（7.14）得：（7.16）

對式(7.15)和(7.16)進行比較系數(shù)后可知,只要（7.17）

成立,格式(7.14)的局部截斷誤差就等于有2階精度第51頁，共93頁，2023年，2月20日，星期三式(7.17)中具有三個未知量,但只有兩個方程,因而有無窮多解。若取,則p=1，這是無窮多解中的一個解，將以上所解的值代入式(7.14)并改寫可得

不難發(fā)現(xiàn)，上面的格式就是改進的歐拉格式。凡滿足條件式（7.17）有一簇形如上式的計算格式，這些格式統(tǒng)稱為二階龍格—庫塔格式。因此改進的歐拉格式是眾多的二階龍格—庫塔法中的一種特殊格式。第52頁，共93頁，2023年，2月20日，星期三若取,則，此時二階龍格-庫塔法的計算公式為

此計算公式稱為變形的二階龍格—庫塔法。式中為區(qū)間的中點。第53頁，共93頁，2023年，2月20日，星期三7.4.3三階龍格-庫塔法

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

為了預(yù)報點的斜率值k3,在區(qū)間內(nèi)有兩個斜率值k1和k2可以用,可將k1,k2加權(quán)平均得出上的平均斜率,從而得到的預(yù)報值第54頁，共93頁，2023年，2月20日，星期三于是可得

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

（7.19）

第55頁，共93頁，2023年，2月20日，星期三7.4.4四階龍格—庫塔法

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

（7.20）

第56頁，共93頁，2023年，2月20日，星期三7.4.5四階龍格—庫塔法算法實現(xiàn)(1)

計算步驟①輸入,h,N②使用龍格—庫塔公式（7.20）計算出y1③輸出，并使轉(zhuǎn)到②直至n>N結(jié)束。第57頁，共93頁，2023年，2月20日，星期三(2)四階龍格—庫塔算法流程圖第58頁，共93頁，2023年，2月20日，星期三程序?qū)崿F(xiàn)(見附錄AA-16四階龍格-庫塔法計算常微分方程初值問題)

例7.4取步長h=0.2，用經(jīng)典格式求解初值問題解:由四階龍格-庫塔公式可得

第59頁，共93頁，2023年，2月20日，星期三可同樣進行其余yi的計算。本例方程的解為，數(shù)值解yi與準確解y(xi)的對照表見教材P163所示，從表中看到所求的數(shù)值解具有4位有效數(shù)字。

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

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

當h值不大時，式中的系數(shù)c可近似地看作為常數(shù)。第61頁，共93頁，2023年，2月20日，星期三然后將步長折半,即以為步長,從節(jié)點xi出發(fā),跨兩步到節(jié)點xi+1,再求得一個近似值,每跨一步的截斷誤差是,因此有這樣由此可得這表明以作為的近似值，其誤差可用步長折半前后兩次計算結(jié)果的偏差來判斷所選步長是否適當?shù)?2頁，共93頁，2023年，2月20日，星期三當要求的數(shù)值精度為ε時：（1）如果Δ>ε，反復(fù)將步長折半進行計算，直至Δ<ε為止,并取其最后一次步長的計算結(jié)果作為（2）如果Δ<ε，反復(fù)將步長加倍，直到Δ>ε為止，并以上一次步長的計算結(jié)果作為。這種通過步長加倍或折半來處理步長的方法稱為變步長法。表面上看，為了選擇步長，每一步都要反復(fù)判斷Δ，增加了計算工作量，但在方程的解y(x)變化劇烈的情況下，總的計算工作量得到減少，結(jié)果還是合算的。第63頁，共93頁，2023年，2月20日，星期三7.5亞當姆斯方法7.5.1亞當姆斯格式龍格-庫塔方法是一類重要算法，但這類算法在每一步都需要先預(yù)報幾個點上的斜率值，計算量比較大?？紤]到計算yi+1之前已得出一系列節(jié)點上的斜率值，能否利用這些已知值來減少計算量呢？這就是亞當姆斯（Adams）方法的設(shè)計思想。

第64頁，共93頁，2023年，2月20日，星期三設(shè)用xi,xi+1兩點的斜率值加權(quán)平均作為區(qū)間上的平均斜率，有計算格式

（7.21）

選取參數(shù)λ，使格式（7.21）具有二階精度。第65頁，共93頁，2023年，2月20日，星期三將在xi處Taylor展開

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

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

第66頁，共93頁，2023年，2月20日，星期三和四階亞當姆斯格式。

(7.22)這里和下面均記上述幾種亞當姆斯格式都是顯式的，算法比較簡單，但用節(jié)點的斜率值來預(yù)報區(qū)間上的平均斜率是個外推過程，效果不夠理想。為了進一步改善精度，變外推為內(nèi)插，即增加節(jié)點xi+1的斜率值來得出上的平均斜率。譬如考察形如第67頁，共93頁，2023年，2月20日，星期三（7.23）

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

可見要使格式(7.23)具有二階精度,需令,這樣就可構(gòu)造二階隱式亞當姆斯格式其實是梯形格式。類似可導(dǎo)出三階隱式亞當姆斯格式第68頁，共93頁，2023年，2月20日，星期三和四階隱式亞當姆斯格式

（7.24）

7.5.2亞當姆斯預(yù)報-校正格式參照改進的歐拉格式的構(gòu)造方法，以四階亞當姆斯為例，將顯式（7.22）和隱式（7.24）相結(jié)合，用顯式公式做預(yù)報，再用隱式公式做校正，可構(gòu)成亞當姆斯預(yù)報-校正格式（7.25）

預(yù)報：

校正：

第69頁，共93頁，2023年，2月20日，星期三這種預(yù)報-校正格式是四步法，它在計算yi+1時不但用到前一步的信息，而且要用到再前面三步的信息，因此它不能自行啟動。在實際計算時，可借助于某種單步法，譬如四階龍格—庫塔法提供開始值。

第70頁，共93頁，2023年，2月20日，星期三例7.5取步長h=0.1,用亞當姆斯預(yù)報-校正公式求解初值問題

的數(shù)值解。解:用四階龍格-庫塔公式求出發(fā)值,計算得：表中的,yi和y(xi)分別為預(yù)報值、校正值和準確解(),以比較計算結(jié)果的精度。再使用亞當姆斯預(yù)報-校正公式(7.25),見教材P166列表算得其余的計算結(jié)果第71頁，共93頁，2023年，2月20日，星期三7.6算法的穩(wěn)定性及收斂性7.6.1穩(wěn)定性穩(wěn)定性在微分方程的數(shù)值解法中是一個非常重要的問題。因為微分方程初值問題的數(shù)值方法是用差分格式進行計算的，而在差分方程的求解過程中，存在著各種計算誤差，這些計算誤差如舍入誤差等引起的擾動，在傳播過程中，可能會大量積累，對計算結(jié)果的準確性將產(chǎn)生影響。這就涉及到算法穩(wěn)定性問題。

第72頁，共93頁，2023年，2月20日，星期三

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

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

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

總體截斷誤差（7.28）

又由于f(x,y)關(guān)于y滿足李普希茨條件,即第77頁，共93頁，2023年，2月20日，星期三代入(7.28)上式，有再利用式（7.27），式（7.28）即上式反復(fù)遞推后，可得（7.29）

第78頁，共93頁，2023年，2月20日，星期三設(shè)（T為常數(shù)）

因為

所以把上式代入式（7.29），得若不計初值誤差，即，則有（7.30）

式(7.30)說明,當時,,即,所以Euler方法是收斂的，且其收斂速度為，即具有一階收斂速度。同時還說明Euler方法的整體截斷誤差為，因此算法的精度為一階。第79頁，共93頁，2023年，2月20日，星期三7.7一階常微分方程組與高階方程我們已介紹了一階常微分方程初值問題的各種數(shù)值解法，對于一階常微分方程組，可類似得到各種解法，而高階常微分方程可轉(zhuǎn)化為一階常微分方程組來求解。

7.7.1一階常微分方程組對于一階常微分方程組的初值問題

（7.31）

可以把單個方程中的f和y看作向量來處理，這樣就可把前面介紹的各種差分算法推廣到求一階方程組初值問題中來。第80頁，共93頁，2023年，2月20日，星期三設(shè)為節(jié)點上的近似解，則有改進的Euler格式為預(yù)報：校正：

（7.32）

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

第81頁，共93頁，2023年，2月20日，星期三式中

（7.34）

第82頁，共93頁，2023年，2月2