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

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

定理1

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

定理2

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

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

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

例1

用歐拉公式求解初值問題

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

xn

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

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

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

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

例2

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

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

xnyn

誤差

xnyn誤差00.20.4

1.1.1840961.34336000.0000880.0017190.60.81.01.4859561.6164761.7378690.0027160.0040240.05818同例1中的歐拉法的計(jì)算結(jié)果比較，明顯改善了精度.第22頁，課件共59頁，創(chuàng)作于2023年2月9.2.4單步法的局部截?cái)嗾`差與階初值問題(1.1),(1.2)的單步法可用一般形式表示為其中多元函數(shù)與f(x,y

)有關(guān)，當(dāng)含有yn+1時(shí)，方法是隱式的，若不含yn+1則為顯式方法，所以顯式單步法可表示為(x,y,h)稱為增量函數(shù)，例如對(duì)歐拉法(2.1)有它的局部截?cái)嗾`差已由(2.3)給出,對(duì)一般顯式單步法則可如下定義.第23頁，課件共59頁，創(chuàng)作于2023年2月

定義1

設(shè)y(x)是初值問題(1.1),(1.2)的準(zhǔn)確解,稱為顯式單步法(2.11)的局部截?cái)嗾`差.

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

所以，局部截?cái)嗾`差可理解為用方法(2.11)計(jì)算一步的誤差，也即公式(2.11)中用準(zhǔn)確解y(x)代替數(shù)值解產(chǎn)生的公式誤差.根據(jù)定義,顯然歐拉法的局部截?cái)嗾`差為第24頁，課件共59頁，創(chuàng)作于2023年2月即為(2.3)的結(jié)果.這里稱為局部截?cái)嗾`差主項(xiàng).顯然Tn+1=O(h2).一般情形的定義如下第25頁，課件共59頁，創(chuàng)作于2023年2月

定義2

設(shè)y(x)是初值問題的準(zhǔn)確解，若存在最大整數(shù)p使顯式單步法(2.11)的局部截?cái)嗾`差滿足則稱方法(2.11)具有p階精度.若將(2.11)展開式寫成則稱為局部截?cái)嗾`差主項(xiàng).第26頁，課件共59頁，創(chuàng)作于2023年2月以上定義對(duì)隱式單步法(2.10)也是適用的.例如，對(duì)后退歐拉法(2.5)其局部截?cái)嗾`差為這里p=1是1階方法，局部截?cái)嗾`差主項(xiàng)為第27頁，課件共59頁，創(chuàng)作于2023年2月同樣對(duì)梯形法(2.7)有所以梯形方法(2.7)是二階的.其局部截?cái)嗾`差主項(xiàng)為第28頁，課件共59頁，創(chuàng)作于2023年2月9.3龍格-庫塔方法

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

第29頁，課件共59頁，創(chuàng)作于2023年2月9.3.1顯式龍格-庫塔法的一般形式上節(jié)給出了顯式單步法的表達(dá)式(2.11),其局部截?cái)嗾`差為O(hp+1)，對(duì)歐拉法Tn+1=O(h2)，即方法為p=1階，若用改進(jìn)的歐拉法(2.9)式，它可表示為此時(shí)增量函數(shù)為第30頁，課件共59頁，創(chuàng)作于2023年2月它比歐拉法的(xn,yn,h)=f(xn,yn),增加了計(jì)算一個(gè)右函數(shù)f的值，可望p=2.若要使得到的公式階數(shù)p更大,就必須包含更多的f

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

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

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

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

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

四階R-K方法的每一步需要計(jì)算四次函數(shù)值f，可以證明其局部截?cái)嗾`差為O(h5).證明從略.第44頁，課件共59頁，創(chuàng)作于2023年2月

例3

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

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

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

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