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

第十三講常微分方程初值問題數(shù)值解法第1頁，共59頁，2023年，2月20日，星期二9.1引言科學(xué)技術(shù)中很多問題都可用微分方程的定解問題來描述，主要有初值問題與邊值問題兩大類，本章只考慮初值問題.常微分方程初值問題中最簡單的例子是人口模型，設(shè)某特定區(qū)域在t0時刻人口為y(t0)=y0已知的，該區(qū)域的人口自然增長率為，人口增長與人口總數(shù)成正比，所以t時刻的人口總數(shù)y(t)滿足以下微分方程第2頁，共59頁，2023年，2月20日，星期二很多物理系統(tǒng)與時間有關(guān)，從衛(wèi)星運行軌道到單擺運動，從化學(xué)反應(yīng)到物種競爭都是隨時間的延續(xù)而不斷變化的.解常微分方程是描述連續(xù)變化的數(shù)學(xué)語言，微分方程的求解就是確定滿足給定方程的可微函數(shù)y(t)，研究它的數(shù)值方法是本章的主要目的.考慮一階常微分方程的初值問題則稱f關(guān)于y滿足利普希茨(Lipschitz)條件，L稱為y的利普希茨常數(shù)(簡稱Lips.常數(shù)).如果存在實數(shù)L>0，使得第3頁，共59頁，2023年，2月20日，星期二

定理1

設(shè)f在區(qū)域D={(x,y)|axb,yR}上連續(xù),關(guān)于y滿足利普希茨條件，則對任意x0[a,b],y0R，常微分方程初值問題(1.1)式和(1.2)式當(dāng)x[a,b]時存在唯一的連續(xù)可微解y(x).解的存在唯一性定理是常微分方程理論的基本內(nèi)容，也是數(shù)值方法的出發(fā)點，此外還要考慮方程的解對擾動的敏感性，它有以下結(jié)論.

定理2

設(shè)f在區(qū)域D(如定理1所定義)

上連續(xù),且關(guān)于y滿足利普希茨條件，設(shè)初值問題的解為y(x,s)，則第4頁，共59頁，2023年，2月20日，星期二這個定理表明解對初值依賴的敏感性，它與右端函數(shù)f有關(guān)，當(dāng)f的Lips.常數(shù)L比較小時，解對初值和右端函數(shù)相對不敏感，可視為好條件.若L較大則可認為壞條件，即病態(tài)問題.如果右端函數(shù)可導(dǎo)，由中值定理有若假定在域D內(nèi)有界,設(shè),則它表明f滿足利普希茨條件，且L的大小反映了右端函數(shù)f關(guān)于y變化的快慢，刻畫了初值問(1.1)式和(1.2)式是否為好條件.這在數(shù)值求解中也是很重要的.第5頁，共59頁，2023年，2月20日，星期二雖然求解常微分方程有各種各樣的解析方法，但解析方法只能用來求解一些特殊類型的方程，實際問題中歸結(jié)出來的微分方程主要靠數(shù)值解法.所謂數(shù)值解法,就是尋求解y(x)在一系列離散節(jié)點上的近似值y1,y2,,yn,yn+1,.相鄰兩個節(jié)點的間距hn=xn+1-xn稱為步長.今后如不特別說明，總是假定hi=h(i=1,2,)為常數(shù),這時節(jié)點為xn=x0+nh(i=0,1,2,)(等距節(jié)點).第6頁，共59頁，2023年，2月20日，星期二初值問題的數(shù)值解法有個基本特點，他們都采取“步進式”，即求解過程順著節(jié)點排列的次序一步一步地向前推進.描述這類算法，只要給出用已知信息yn,yn-1,yn-2,計算yn+1的遞推公式.本章首先要對常微分方程(1.1)離散化，建立求解數(shù)值解的遞推公式.一類是計算yn+1時只用到前一點的值yn，稱為單步法.另一類是用到y(tǒng)n+1前面k點的值yn,yn-1,,yn-k+1，稱為k步法.其次，要研究公式的局部截斷誤差和階，數(shù)值解yn與精確解y(xn)的誤差估計及收斂性，還有遞推公式的計算穩(wěn)定性等問題.第7頁，共59頁，2023年，2月20日，星期二9.2簡單的數(shù)值方法9.2.1歐拉法與后退歐拉法我們知道，在xy平面上，微分方程(1.1)式的解y=f(x)稱作它的積分曲線，積分曲線上一點(x,y)的切線斜率等于函數(shù)f(x,y)的值.如果按f(x,y)在xy平面上建立一個方向場，那么，積分曲線上每一點的切線方向均與方向場在該點的方向相一致.基于上述幾何解釋，我們從初始點P0(x0,y0)出發(fā),先依方向場在該點的方向推進到x=x1上一點P1，然后再從P1點依方向場在該點的方向推進到x=x2

上一點P2,循環(huán)前進做出一條折線P0P1P2.第8頁，共59頁，2023年，2月20日，星期二一般地，設(shè)已做出該折線的頂點Pn,過Pn(xn,yn)依方向場的方向再推進到Pn+1(xn+1,yn+1)，顯然兩個頂點Pn,Pn+1的坐標有關(guān)系這就是著名的(顯式)歐拉(Euler)公式.若初值y0已知，則依公式(2.1)可逐次逐步算出各點數(shù)值解.即第9頁，共59頁，2023年，2月20日，星期二

例1

用歐拉公式求解初值問題

解取步長h=0.1，歐拉公式的具體形式為其中xn=nh=0.1n(n=0,1,,10),已知y0=1,由此式可得第10頁，共59頁，2023年，2月20日，星期二依次計算下去，部分計算結(jié)果見下表.與準確解相比，可看出歐拉公式的計算結(jié)果精度很差.

xn

歐拉公式數(shù)值解yn準確解y(xn)

誤差

0.20.40.60.81.0

1.1918181.3582131.5089661.6497831.784770

1.1832161.3416411.4832401.6124521.732051

0.0086020.0165720.0257260.0373310.052719第11頁，共59頁，2023年，2月20日，星期二歐拉公式具有明顯的幾何意義,就是用折線近似代替方程的解曲線，因而常稱公式(2.1)為歐拉折線法.還可以通過幾何直觀來考察歐拉方法的精度.假設(shè)yn=y(xn),即頂點Pn落在積分曲線y=y(x)上，那么，按歐拉方法做出的折線PnPn+1便是y=y(x)過點Pn的切線.從圖形上看,這樣定出的頂點Pn+1顯著地偏離了原來的積分曲線，可見歐拉方法是相當(dāng)粗糙的.第12頁，共59頁，2023年，2月20日，星期二為了分析計算公式的精度，通?？捎锰├照归_將y(xn+1)在xn處展開，則有在yn=y(xn)的前提下，f(xn,yn)=f(xn,y(xn))=y(xn).于是可得歐拉法(2.1)的公式誤差為稱為此方法的局部截斷誤差.第13頁，共59頁，2023年，2月20日，星期二如果對方程(1.1)從xn到xn+1積分，得右端積分用左矩形公式hf(xn,y(xn))近似，再以yn代替y(xn)，yn+1代替y(xn+1)也得到歐拉公式(2.1)，局部截斷誤差也是(2.3).稱為(隱式)后退的歐拉公式.它也可以通過利用均差近似導(dǎo)數(shù)y(xn+1)，即如果右端積分用右矩形公式hf(xn+1,y(xn+1))近似，則得到另一個公式直接得到.第14頁，共59頁，2023年，2月20日，星期二后退的歐拉公式與歐拉公式有著本質(zhì)的區(qū)別,后者是關(guān)于yn+1的一個直接計算公式，這類公式稱作是顯式的；前者公式的右端含有未知的yn+1，它實際上是關(guān)于yn+1的一個函數(shù)方程,這類方程稱作是隱式的.顯式與隱式兩類方法各有特點，考慮到數(shù)值穩(wěn)定性等其他因素，人們有時需要選用隱式方法，但使用顯式算法遠比隱式方便.隱式方程(2.5)通常用迭代法求解，而迭代過程的實質(zhì)是逐步顯式化.第15頁，共59頁，2023年，2月20日，星期二設(shè)用歐拉公式給出迭代初值

，用它代入(2.5)式的右端，使之轉(zhuǎn)化為顯式，直接計算得然后再用代入(2.5)式，又有如此反復(fù)進行，得第16頁，共59頁，2023年，2月20日，星期二由于f(x,y)對y滿足Lipschitz條件(1.3).由(2.6)減(2.5)得由此可知，只要hL<1，迭代法(2.6)就收斂到解yn+1.關(guān)于后退歐拉方法的公式誤差，從積分公式看到它與歐拉法是相似的.第17頁，共59頁，2023年，2月20日，星期二9.2.2梯形方法為得到比歐拉法精度高的計算公式，在等式(2.4)

右端積分用梯形求積公式近似,并用yn代替y(xn),yn+1代替y(xn+1)，則得稱為矩形方法.矩形方法是隱式單步法，用迭代法求解，同后退的歐拉方法一樣，仍用歐拉法提供迭代初值，則矩形迭代公式為第18頁，共59頁，2023年，2月20日，星期二為了分析迭代過程的收斂性,將(2.7)與(2.8)相減,得于是有使得則當(dāng)k→∞時有,這說明迭代過程(2.8)是收斂的.第19頁，共59頁，2023年，2月20日，星期二9.2.3改進的歐拉公式我們看到，梯形方法雖然提高了精度，但其算法復(fù)雜，在應(yīng)用迭代公式(2.8)進行實際計算時，每迭代一次，都要重新計算函數(shù)f(x,y

)的值，而迭代又要反復(fù)進行若干次，計算量很大，而且往往難以預(yù)測.為了控制計算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計算，這就簡化了算法.具體地說，我們先用歐拉公式求得一個初步的近似值，稱之為預(yù)測值，此預(yù)測值的精度可能很差，再用梯形公式(2.7)將它校正一次，即按(2.8)式迭代一次得yn+1，這個結(jié)果稱之為校正值.第20頁，共59頁，2023年，2月20日，星期二這樣建立的預(yù)測-校正系統(tǒng)通常稱為改進的歐拉公式:或表示為下列平均化形式(2.9)預(yù)測校正第21頁，共59頁，2023年，2月20日，星期二

例2

用改進的歐拉法解例1中的初值問題(2.2).

解仍取步長h=0.1,改進的歐拉公式為部分計算結(jié)果見下表

xnyn

誤差

xnyn誤差00.20.4

1.1.1840961.34336000.0000880.0017190.60.81.01.4859561.6164761.7378690.0027160.0040240.05818同例1中的歐拉法的計算結(jié)果比較，明顯改善了精度.第22頁，共59頁，2023年，2月20日，星期二9.2.4單步法的局部截斷誤差與階初值問題(1.1),(1.2)的單步法可用一般形式表示為其中多元函數(shù)與f(x,y

)有關(guān)，當(dāng)含有yn+1時，方法是隱式的，若不含yn+1則為顯式方法，所以顯式單步法可表示為(x,y,h)稱為增量函數(shù)，例如對歐拉法(2.1)有它的局部截斷誤差已由(2.3)給出,對一般顯式單步法則可如下定義.第23頁，共59頁，2023年，2月20日，星期二

定義1

設(shè)y(x)是初值問題(1.1),(1.2)的準確解,稱為顯式單步法(2.11)的局部截斷誤差.

Tn+1之所以稱為局部的，是假設(shè)在xn前各步?jīng)]有誤差.當(dāng)yn=y(xn)時，計算一步，則有

所以，局部截斷誤差可理解為用方法(2.11)計算一步的誤差，也即公式(2.11)中用準確解y(x)代替數(shù)值解產(chǎn)生的公式誤差.根據(jù)定義,顯然歐拉法的局部截斷誤差為第24頁，共59頁，2023年，2月20日，星期二即為(2.3)的結(jié)果.這里稱為局部截斷誤差主項.顯然Tn+1=O(h2).一般情形的定義如下第25頁，共59頁，2023年，2月20日，星期二

定義2

設(shè)y(x)是初值問題的準確解，若存在最大整數(shù)p使顯式單步法(2.11)的局部截斷誤差滿足則稱方法(2.11)具有p階精度.若將(2.11)展開式寫成則稱為局部截斷誤差主項.第26頁，共59頁，2023年，2月20日，星期二以上定義對隱式單步法(2.10)也是適用的.例如，對后退歐拉法(2.5)其局部截斷誤差為這里p=1是1階方法，局部截斷誤差主項為第27頁，共59頁，2023年，2月20日，星期二同樣對梯形法(2.7)有所以梯形方法(2.7)是二階的.其局部截斷誤差主項為第28頁，共59頁，2023年，2月20日，星期二9.3龍格-庫塔方法

對許多實際問題來說，歐拉公式與改進歐拉公式精度還不能滿足要求，為此從另一個角度來分析這兩個公式的特點，從而探索一條構(gòu)造高精度方法的途徑.

第29頁，共59頁，2023年，2月20日，星期二9.3.1顯式龍格-庫塔法的一般形式上節(jié)給出了顯式單步法的表達式(2.11),其局部截斷誤差為O(hp+1)，對歐拉法Tn+1=O(h2)，即方法為p=1階，若用改進的歐拉法(2.9)式，它可表示為此時增量函數(shù)為第30頁，共59頁，2023年，2月20日，星期二它比歐拉法的(xn,yn,h)=f(xn,yn),增加了計算一個右函數(shù)f的值，可望p=2.若要使得到的公式階數(shù)p更大,就必須包含更多的f

值.實際上從方程(1.1)等價的積分形式(2.4)

，即若要使公式階數(shù)提高，就必須使右端積分的數(shù)值求積公式精度提高，它必然要增加求積節(jié)點，為此可將(3.3)的右端用求積公式表示為第31頁，共59頁，2023年，2月20日，星期二一般說來，點數(shù)r越多，精度越高，上式右端相當(dāng)于增量函數(shù)(x,y,h)，為得到便于計算的顯式方法,可類似于改進歐拉法(3.1),(3.2)，將公式表示為其中這里ci,i,ij均為常數(shù).(3.4)和(3.5)稱為r級顯式龍格-庫塔(Runge-Kutta)法,簡稱R-K方法.第32頁，共59頁，2023年，2月20日，星期二當(dāng)r=1,(xn,yn,h)=f(xn,yn)時，就是歐拉法，此時方法的階為p=1.當(dāng)r=2時，改進歐拉法(3.1)式就是其中的一種,下面將證明其階p=2.要使公式(3.4),(3.5)具有更高的階p，就要增加點數(shù)r.下面我們只就r=2推導(dǎo)R-K方法.并給出r=3,4時的常用公式，其推導(dǎo)方法與r=2時類似，只是推導(dǎo)計算較復(fù)雜.第33頁，共59頁，2023年，2月20日，星期二9.3.2二階顯式R-K方法對r=2的R-K方法，由(3.4),(3.5)式可得如下計算公式這里c1,c2,2,21均為待定常數(shù)，我們希望適當(dāng)選取這些系數(shù)，使公式階數(shù)p盡量高.根據(jù)局部截斷誤差定義，推導(dǎo)出(3.6)的局部截斷誤差為第34頁，共59頁，2023年，2月20日，星期二其中這里yn=y(xn),yn+1=y(xn+1).為得到Tn+1的階p，要將上式各項在(xn,

yn)處做泰勒展開，由于f(x,y

)是二元函數(shù)，故要用二元泰勒展開，各項展開式為第35頁，共59頁，2023年，2月20日，星期二將以上結(jié)果代入(3.7)，則有第36頁，共59頁，2023年，2月20日，星期二要使公式(3.6)具有p=2階，必須使即(3.9)的解是不唯一的.可令c2=a≠0，則得這樣得到的公式稱為二階R-K方法.第37頁，共59頁，2023年，2月20日，星期二則由此可以看出在改進的歐拉公式中相當(dāng)于取(xn,yn),(xn+1,yn+1)兩點處斜率的平均值，近似代替平均斜率，其精度比歐拉公式提高了.如取a=1/2，則c1=c2=1/2,2=21=1.這就是改進的歐拉公式(3.1)式.第38頁，共59頁，2023年，2月20日，星期二稱為中點公式(變形的歐拉公式)，相當(dāng)于數(shù)值積分的中矩形公式.也可以表示為如取a=1，則c1=0,c2=1,2=21=1/2.得計算公式第39頁，共59頁，2023年，2月20日，星期二對r=2的R-K公式(3.6)能否使局部誤差提高到O(h4)?為此需把K2多展開一項，從(3.8)的看到展開式中的項是不能通過選擇參數(shù)消掉的，實際上要使h3

的項為零，需增加3個方程，要確定4個參數(shù)c1,c2,2及21，這是不可能的.故r=2的顯式R-K方法的階只能是p=2，而不能得到三階公式.第40頁，共59頁，2023年，2月20日，星期二9.3.3三階與四階顯式R-K方法要得到三階顯式R-K方法，必須r=3.此時計算(3.4),(3.5)的公式表示為其中c1,c2,c3及2,21,3,31,32均為待定常數(shù)，公式(3.11)的局部截斷誤差為第41頁，共59頁，2023年，2月20日，星期二只要K1,K2將按二元泰勒展開，使Tn+1＝O(h4)，可得待定參數(shù)滿足方程第42頁，共59頁，2023年，2月20日，星期二這是8個未知數(shù)6個方程的方程組，解不是唯一的.可以得到很多公式.滿足條件(3.12)的公式(3.11)統(tǒng)稱為三階R-K公式.下面只給出其中一個常見的公式.此公式稱為庫塔三階方法.第43頁，共59頁，2023年，2月20日，星期二繼續(xù)上述過程,經(jīng)過較復(fù)雜的數(shù)學(xué)演算,可以導(dǎo)出各種四階R-K公式,下列經(jīng)典公式是其中常用的一個：

四階R-K方法的每一步需要計算四次函數(shù)值f，可以證明其局部截斷誤差為O(h5).證明從略.第44頁，共59頁，2023年，2月20日，星期二

例3

設(shè)取步長h=0.2,從x=0直到x=1用四階龍格-庫塔方法求解例1中的初值問題(2.2)式.

解這里,經(jīng)典的四階龍格-庫塔公式(3.13)具有第45頁，共59頁，2023年，2月20日，星期二計算結(jié)果見下表

xnyny(xn)

xnyny(xn)

00.20.4

1.1.18321.3417

01.18323.24160.60.81.0

1.48331.61251.7321

1.48321.61251.7321比較例3和例2的計算結(jié)果，顯然以龍格-庫塔方法的精度為高.要注意，雖然四階龍格-庫塔方法的計算量(每一步要4次計算函數(shù)f)比改進的歐拉方法(它是一種二階龍格-庫塔方法，每一步只要2次計算函數(shù)f)大一倍，但由于這里放大了步長(h=0.2)，兩個例子的算法所耗費的計算量幾乎相同.這個例子又一次顯示了選擇算法的重要意義.第46頁，共59頁，2023年，2月20日，星期二然而值得指出的是，龍格-庫塔方法的推導(dǎo)基于泰勒展開方法，因而它要求所求的解具有較好的光滑性質(zhì).反之，如果解的光滑性差，那么，使用龍格-庫塔方法求得的數(shù)值解，其精度可能反而不如改進的歐拉方法.實際計算時，我們應(yīng)當(dāng)針對問題的具體特點選擇合適的算法.第47頁，共59頁，2023年，2月20日，星期二微分方程組和高階方程初值問題的數(shù)值解歐拉方法和龍格-庫塔方法可直接推廣到微分方程組向前歐拉公式高階方程需要先降階為一階微分方程組

第48頁，共59頁，2023年，2月20日，星期二龍格—庫塔方法的MATLAB實現(xiàn)[t,x]=ode23(@f,ts,x0,opt)3級2階龍格-庫塔公式[t,x]=ode45(@f,ts,x0,opt)5級4階龍格-庫塔公式f是待解方程寫成的函數(shù)m文件：functiondx=f(t,x)dx=[f1;f2;;fn];ts=[t0,t1,…,tf]輸出指定時刻t0,t1,…,tf的函數(shù)值ts=t0:k:tf輸出[t0,tf]內(nèi)等分點處的函數(shù)值x0為函數(shù)初值(n維)輸出t=ts,x為相應(yīng)函數(shù)值(n維)opt為選項，缺省時精度為：相對誤差10-3，絕對誤差10-6,計算步長按精度要求自動調(diào)整第49頁，共59頁，2023年，2月20日，星期二微分方程（組）的數(shù)值解

事實上，能夠求得解析解的微分方程或微分方程組少之又少，多數(shù)情況下需要求出微分方程（組）的數(shù)值解。Matlab中求微分方程數(shù)值解的函數(shù)有五個：ode45，ode23，ode113，ode15s，ode23s。調(diào)用格式為[t,x]=solver(‘f’,ts,x0,options)第50頁，共59頁，2023年，2月20日，星期二需要特別注意的是：①solver可以取以上五個函數(shù)之一，不同的函數(shù)代表不同的內(nèi)部算法：ode23運用組合的2/3階龍格—庫塔—費爾貝算法，ode45運用組合的4/5階龍格—庫塔—費爾貝算法。通常使用函數(shù)ode45；②f是由待解方程寫成的m文件的文件名；③ts=[t0,tf]，t0、tf為自變量的初值和終值；④x0為函數(shù)的初值；第51頁，共59頁，2023年，2月20日，星期二⑤options用于設(shè)