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

第九章常微分方程初值問(wèn)題數(shù)值解法

數(shù)值分析2022/10/201NumericalAnalysis第九章常微分方程初值問(wèn)題數(shù)值解法

數(shù)值分析2022/10本章內(nèi)容簡(jiǎn)單的數(shù)值方法歐拉法與后退歐拉法梯形方法改進(jìn)歐拉公式單步法的局部截?cái)嗾`差與階龍格-庫(kù)塔方法顯式龍格-庫(kù)塔法的一般形式二階顯式R-K方法三階與四階顯式R-K方法2022/10/202NumericalAnalysis本章內(nèi)容簡(jiǎn)單的數(shù)值方法2022/10/152Numerica9.1引言科學(xué)技術(shù)中常常需要求解常微分方程的定解問(wèn)題.這類(lèi)問(wèn)題最簡(jiǎn)單的形式，是本章將著重考察的一階方程的初值問(wèn)題我們知道，只有f(x,y)適當(dāng)光滑—譬如關(guān)于y滿足利普希茨(Lipschitz)條件理論上就可以保證初值問(wèn)題的解y＝f(x)存在并且唯一.2022/10/203NumericalAnalysis9.1引言科學(xué)技術(shù)中常常需要雖然求解常微分方程有各種各樣的解析方法，但解析方法只能用來(lái)求解一些特殊類(lèi)型的方程，實(shí)際問(wèn)題中歸結(jié)出來(lái)的微分方程主要靠數(shù)值解法.所謂數(shù)值解法,就是尋求解y(x)在一系列離散節(jié)點(diǎn)上的近似值y1,y2,,yn,yn+1,.相鄰兩個(gè)節(jié)點(diǎn)的間距hn=xn+1-xn稱(chēng)為步長(zhǎng).今后如不特別說(shuō)明，總是假定hi=h(i=1,2,)為定數(shù),這時(shí)節(jié)點(diǎn)為xn=x0+nh(i=0,1,2,)(等距節(jié)點(diǎn)).2022/10/204NumericalAnalysis雖然求解常微分方程有各種各樣的解析方法，但解初值問(wèn)題的數(shù)值解法有個(gè)基本特點(diǎn)，他們都采取“步進(jìn)式”，即求解過(guò)程順著節(jié)點(diǎn)排列的次序一步一步地向前推進(jìn).描述這類(lèi)算法，只要給出用已知信息yn,yn-1,yn-2,計(jì)算yn+1的遞推公式.首先，要對(duì)微分方程離散化，建立求解數(shù)值解的遞推公式.一類(lèi)是計(jì)算yn+1時(shí)只用到前一點(diǎn)的值yn，稱(chēng)為單步法.另一類(lèi)是用到y(tǒng)n+1前面k點(diǎn)的值yn,yn-1,,yn-k+1，稱(chēng)為k步法.其次，要研究公式的局部截?cái)嗾`差和階，數(shù)值解yn與精確解y(xn)的誤差估計(jì)及收斂性，還有遞推公式的計(jì)算穩(wěn)定性等問(wèn)題.2022/10/205NumericalAnalysis初值問(wèn)題的數(shù)值解法有個(gè)基本特點(diǎn)，他們都采取9.2簡(jiǎn)單的數(shù)值方法與基本概念9.2.1歐拉法與后退歐拉法我們知道，在xy平面上，微分方程(1.1)式的解y=f(x)稱(chēng)作它的積分曲線，積分曲線上一點(diǎn)(x,y)的切線斜率等于函數(shù)f(x,y)的值.如果按f(x,y)在xy平面上建立一個(gè)方向場(chǎng)，那么，積分曲線上每一點(diǎn)的切線方向均與方向場(chǎng)在該點(diǎn)的方向相一致.基于上述幾何解釋?zhuān)覀儚某跏键c(diǎn)P0(x0,y0)出發(fā),先依方向場(chǎng)在該點(diǎn)的方向推進(jìn)到x=x1上一點(diǎn)P1，然后再?gòu)腜1點(diǎn)依方向場(chǎng)在該點(diǎn)的方向推進(jìn)到x=x2

上一點(diǎn)P2,循環(huán)前進(jìn)做出一條折線P0P1P2.2022/10/206NumericalAnalysis9.2簡(jiǎn)單的數(shù)值方法與基本概念9.2.1歐拉法與后一般地，設(shè)已做出該折線的頂點(diǎn)Pn,過(guò)Pn(xn,yn)依方向場(chǎng)的方向再推進(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ù)值解.即2022/10/207NumericalAnalysis一般地，設(shè)已做出該折線的頂點(diǎn)Pn,過(guò)Pn(xn,y

例1用歐拉公式求解初值問(wèn)題

解取步長(zhǎng)h=0.1，歐拉公式的具體形式為其中xn=nh=0.1n(n=0,1,,10),已知y0=1,由此式可得2022/10/208NumericalAnalysis例1用歐拉公式求解初值問(wèn)題解取步長(zhǎng)h=0依次計(jì)算下去，部分計(jì)算結(jié)果見(jiàn)下表.與準(zhǔn)確解相比，可看出歐拉公式的計(jì)算結(jié)果精度很差.

xn

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

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

，用它代入(2.5)式的右端，使之轉(zhuǎn)化為顯式，直接計(jì)算得然后再用代入(2.5)式，又有如此反復(fù)進(jìn)行，得2022/10/2014NumericalAnalysis設(shè)用歐拉公式給出迭代初值，用它代入(2.5)由于f(x,y)對(duì)y滿足Lipschitz條件(1.3).由(2.6)減(2.5)得由此可知，只要hL<1，迭代法(2.6)就收斂到解.關(guān)于后退歐拉方法的公式誤差，從積分公式看到它與歐拉法是相似的.2022/10/2015NumericalAnalysis由于f(x,y)對(duì)y滿足Lipschitz條件(1.3).9.2.2梯形方法為得到比歐拉法精度高的計(jì)算公式，在等式(2.4)右端積分用梯形求積公式近似,并用yn代替y(xn),yn+1代替y(xn+1)，則得稱(chēng)為矩形方法.矩形方法是隱式單步法，用迭代法求解，同后退的歐拉方法一樣，仍用歐拉法提供迭代初值，則矩形迭代公式為2022/10/2016NumericalAnalysis9.2.2梯形方法為得到比歐拉法精度高的計(jì)算公為了分析迭代過(guò)程的收斂性,將(2.7)與(2.8)相減,得于是有使得則當(dāng)k→∞時(shí)有,這說(shuō)明迭代過(guò)程(2.8)是收斂的.2022/10/2017NumericalAnalysis為了分析迭代過(guò)程的收斂性,將(2.7)與(2.8)相減,9.2.3單步法的局部截?cái)嗾`差與階初值問(wèn)題(1.1),(1.2)的單步法可用一般形式表示為其中多元函數(shù)與f(x,y

)有關(guān)，當(dāng)含有yn+1時(shí)，方法是隱式的，若不含yn+1則為顯式方法，所以顯式單步法可表示為(x,y,h)稱(chēng)為增量函數(shù)，例如對(duì)歐拉法(2.1)有它的局部截?cái)嗾`差(2.3)已由給出,對(duì)一般顯式單步法則可如下定義.2022/10/2018NumericalAnalysis9.2.3單步法的局部截?cái)嗾`差與階初值問(wèn)題(1

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

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

所以，局部截?cái)嗾`差可理解為用方法(2.10)計(jì)算一步的誤差，也即公式(2.10)中用準(zhǔn)確解y(x)代替數(shù)值解產(chǎn)生的公式誤差.根據(jù)定義,顯然歐拉法的局部截?cái)嗾`差為2022/10/2019NumericalAnalysis定義1設(shè)y(x)是初值問(wèn)題(1.1),(1.2)的即為(2.3)的結(jié)果.這里稱(chēng)為局部截?cái)嗾`差主項(xiàng).顯然Tn+1=O(h2).一般情形的定義如下

定義2設(shè)y(x)是初值問(wèn)題的準(zhǔn)確解，若存在最大整數(shù)p使顯式單步法(2.10)的局部截?cái)嗾`差滿足則稱(chēng)方法(2.10)具有p階精度.2022/10/2020NumericalAnalysis即為(2.3)的結(jié)果.這里稱(chēng)為局部截?cái)嗾`若將(2.10)展開(kāi)式寫(xiě)成則稱(chēng)為局部截?cái)嗾`差主項(xiàng).以上定義對(duì)隱式單步法(2.9)也是適用的.例如，對(duì)后退歐拉法(2.5)其局部截?cái)嗾`差為這里p=1是1階方法，局部截?cái)嗾`差主項(xiàng)為2022/10/2021NumericalAnalysis若將(2.10)展開(kāi)式寫(xiě)成則同樣對(duì)梯形法(2.7)有所以梯形方法(2.7)是二階的.其局部截?cái)嗾`差主項(xiàng)為2022/10/2022NumericalAnalysis同樣對(duì)梯形法(2.7)有所以梯形方法(2.7)是二階的.9.2.4改進(jìn)的歐拉公式我們看到，梯形方法雖然提高了精度，但其算法復(fù)雜，在應(yīng)用迭代公式(2.9)進(jìn)行實(shí)際計(jì)算時(shí)，每迭代一次，都要重新計(jì)算函數(shù)f(x,y

)的值，而迭代又要反復(fù)進(jìn)行若干次，計(jì)算量很大，而且往往難以預(yù)測(cè).為了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計(jì)算，這就簡(jiǎn)化了算法.具體地說(shuō)，我們先用歐拉公式求得一個(gè)初步的近似值，稱(chēng)之為預(yù)測(cè)值，此預(yù)測(cè)值的精度可能很差，再用梯形公式(2.7)將它校正一次，即按(2.8)式迭代一次，這個(gè)結(jié)果稱(chēng)之為校正值.2022/10/2023NumericalAnalysis9.2.4改進(jìn)的歐拉公式我們看到，梯形方法雖然提高這樣建立的預(yù)測(cè)—校正系統(tǒng)通常稱(chēng)為改進(jìn)的歐拉公式:或表為下列平均化形式(2.13)預(yù)測(cè)校正2022/10/2024NumericalAnalysis這樣建立的預(yù)測(cè)—校正系統(tǒng)通常稱(chēng)為改進(jìn)的歐拉公式:或表為下列平

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

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

xnyn

誤差

xnyn誤差00.20.41.1.1840961.34336000.0000880.0017190.60.81.01.4859561.6164761.7378690.0027160.0040240.05818同例1中的歐拉法的計(jì)算結(jié)果比較，明顯改善了精度.2022/10/2025NumericalAnalysis例2用改進(jìn)的歐拉法解例1中的初值問(wèn)題(2.2).

例

(兩種方法的精度比較)用歐拉公式和改進(jìn)的歐拉公式解下述初值問(wèn)題，并與其準(zhǔn)確解y=e-x+x進(jìn)行比較.解取步長(zhǎng)h=0.1，xk=kh(k=0,1,,6).用兩種方法進(jìn)行計(jì)算對(duì)應(yīng)結(jié)果及絕對(duì)誤差見(jiàn)下表2022/10/2026NumericalAnalysis例(兩種方法的精度比較)用歐拉公式和改進(jìn)的

xn歐拉公式改進(jìn)歐拉公式

yn

誤差

yn

誤差0.00.10.20.20.40.50.61.0000001.0000001.0100001.0290001.0561001.0904901.13144104.8×10-38.7×10-31.2×10-31.4×10-31.6×10-31.8×10-311.0050001.0190251.0412181.0708021.1070761.14940401.6×10-32.9×10-34.0×10-34.8×10-35.5×10-35.9×10-32022/10/2027NumericalAnalysis歐拉公式改進(jìn)歐拉公式y(tǒng)n誤差y9.3龍格—庫(kù)塔方法

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

2022/10/2028NumericalAnalysis9.3龍格—庫(kù)塔方法對(duì)許多實(shí)際問(wèn)題來(lái)說(shuō)9.3.1顯式龍格—庫(kù)塔法的一般形式上節(jié)給出了顯式單步法的表達(dá)式(2.10),其局部截?cái)嗾`差為O(hp+1)，對(duì)歐拉法Tn+1=O(h2)，即方法為p=1階，若用改進(jìn)歐拉法(2.13)，它可表為此時(shí)增量函數(shù)為2022/10/2029NumericalAnalysis9.3.1顯式龍格—庫(kù)塔法的一般形式上節(jié)給出了它比歐拉法的(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)的右端用求積公式表示為2022/10/2030NumericalAnalysis它比歐拉法的(xn,yn,h)=f(xn,yn),一般說(shuō)來(lái)，點(diǎn)數(shù)r越多，精度越高，上式右端相當(dāng)于增量函數(shù)(xn,yn,h)，為得到便于計(jì)算的顯式方法，可類(lèi)似于改進(jìn)歐拉法(3.1),(3.2)，將公式表示為其中這里ci,λi,μij均為常數(shù).(3.4)和(3.5)稱(chēng)為r級(jí)顯式龍格-庫(kù)塔(Runge-Kutta)法,簡(jiǎn)稱(chēng)R-K方法.2022/10/2031NumericalAnalysis一般說(shuō)來(lái)，點(diǎn)數(shù)r越多，精度越高，上式右端相當(dāng)于增量函數(shù)(當(dāng)r=1,(xn,yn,h)=f(xn,yn)時(shí)，就是歐拉法，此時(shí)方法的階為p=1.當(dāng)r=2時(shí)，改進(jìn)歐拉法(3.1),(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í)類(lèi)似，只是計(jì)算較復(fù)雜.2022/10/2032NumericalAnalysis當(dāng)r=1,(xn,yn,h)=f(x9.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)嗾`差為2022/10/2033NumericalAnalysis9.3.2二階顯式R-K方法對(duì)r=2其中這里yn=y(xn),yn+1=y(xn+1).為得到Tn+1的階p，要將上式各項(xiàng)在(xn,

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

)是二元函數(shù)，故要用二元泰勒展開(kāi)，各項(xiàng)展開(kāi)式為2022/10/2034NumericalAnalysis其中這里yn=y(xn),yn+1=y(xn+1).為將以上結(jié)果代入(3.7)，則有2022/10/2035NumericalAnalysis將以上結(jié)果代入(3.7)，則有2022/10/1535Num要使公式(3.6)具有p=2階，必須使即(3.9)的解是不唯一的.可令c2=a≠0，則得這樣得到的公式稱(chēng)為二階R-K方法.2022/10/2036NumericalAnalysis要使公式(3.6)具有p=2階，必須使即(3.9)的解是不唯則由此可以看出在改進(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).2022/10/2037NumericalAnalysis則由此可以看出在改進(jìn)的歐拉公式中相當(dāng)于取(xn,y稱(chēng)為中點(diǎn)公式(變形的歐拉公式)，相當(dāng)于數(shù)值積分的中矩形公式.也可以表示為如取a=1，則c1=0,c2=1,λ2=μ21=1/2.得計(jì)算公式2022/10/2038NumericalAnalysis稱(chēng)為中點(diǎn)公式(變形的歐拉公式)，相當(dāng)于數(shù)值積分的中矩形公式.對(duì)r=2的R-K公式(3.6)能否使局部誤差提高到O(h4)?為此需把K2多展開(kāi)一項(xiàng)，從(3.8)的看到展開(kāi)式中的項(xiàng)是不能通過(guò)選擇參數(shù)削掉的，實(shí)際上要使h3

的項(xiàng)為零，需增加3個(gè)方程，要確定4個(gè)參數(shù)c1,c2,λ2,μ21，這是不可能的.故r＝2的顯式R-K方法的階只能是p=2，而不能得到三階公式.2022/10/2039NumericalAnalysis對(duì)r=2的R-K公式(3.6)能否使局部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)嗾`差為2022/10/2040NumericalAnalysis9.3.3三階與四階顯式R-K方法要只要K1,K2將按二元泰勒展開(kāi)，使Tn+1＝O(h4)，可得待定參數(shù)滿足方程2022/10/2041NumericalAnalysis只要K1,K2將按二元泰勒展開(kāi)，使Tn+1＝O(h4)，可這是8個(gè)未知數(shù)6個(gè)方程的方程組，解不是唯一的.可以得到很多公式.滿足條件(3.12)的公式(3.11)統(tǒng)稱(chēng)為三階R-K公式.下面只給出其中一個(gè)常見(jiàn)的公式.此公式稱(chēng)為庫(kù)塔三階方法.2022/10/2042NumericalAnalysis這是8個(gè)未知數(shù)6個(gè)方程的方程組，解不是唯一的.可以得到很多繼續(xù)上述過(guò)程，經(jīng)過(guò)較復(fù)雜的數(shù)學(xué)演算，可以導(dǎo)出各種四階R-K公式，下列經(jīng)典公式是其中常用的一個(gè)：

四階R-K方法的每一步需要計(jì)算四次函數(shù)值f，可以證明其局部截?cái)嗾`差為O(h5).(例3見(jiàn)書(shū)p349)2022/10/2043NumericalAnalysis繼續(xù)上述過(guò)程，經(jīng)過(guò)較復(fù)雜的數(shù)學(xué)演算，可以導(dǎo)出然而值得指出的是，龍格-庫(kù)塔方法的推導(dǎo)基于泰勒展開(kāi)方法，因而它要求所求的解具有較好的光滑性質(zhì).反之，如果解的光滑性差，那么，使用龍格-庫(kù)塔方法求得的數(shù)值解，其精度可能反而不如改