數(shù)值分析第九章常微分方程數(shù)值解法

1第九章常微分方程數(shù)值解法§1Euler方法§2改進(jìn)的Euler方法§3Runge-Kutta法§4線(xiàn)性多步法2近似求解3一階常微分方程初值問(wèn)題為保證上述問(wèn)題解的存在唯一性,對(duì)函數(shù)f

做如下假設(shè):f關(guān)于y的偏導(dǎo)數(shù)有界,即存在常數(shù)L>0,使得本章中的L均指此數(shù).常微分方程解的存在唯一性4數(shù)值解法:

就是尋求解y(x)在一系列離散點(diǎn)上的近似值相鄰兩個(gè)節(jié)點(diǎn)的間距稱(chēng)為步長(zhǎng).本章均假設(shè)此時(shí)節(jié)點(diǎn)常微分方程數(shù)值解法5

初值問(wèn)題的各種差分方法都采用“步進(jìn)式”，即求解過(guò)程順著節(jié)點(diǎn)排列的次序一步一步地向前推進(jìn).只要給出從已知信息計(jì)算的遞推公式，此類(lèi)計(jì)算格式統(tǒng)稱(chēng)為差分格式.常微分方程差分格式6數(shù)值求解一階常微分方程初值問(wèn)題難點(diǎn):如何離散y

?常見(jiàn)離散方法差商近似導(dǎo)數(shù)數(shù)值積分方法

Taylor展開(kāi)方法7結(jié)果有設(shè)用y(xn)的近似值yn代入上式右端，記所求結(jié)果為yn+1，這樣導(dǎo)出的計(jì)算公式將方程y=f(x,y)中的y在xn點(diǎn)用向前差商近似差商近似導(dǎo)數(shù)(Euler公式)Euler公式,單步法,顯式公式8Euler公式計(jì)算總結(jié)設(shè)的精確解為將區(qū)間[a,b]N等分,分點(diǎn)xn=x0+nh(n=1,2,…,N)步長(zhǎng)用Euler公式求得yn(n=1,2,…,N.)9結(jié)果有由此得差分格式計(jì)算公式將方程y=f(x,y)中的y在xn+1點(diǎn)用向后差商近似向后Euler公式,單步法,隱式公式差商近似導(dǎo)數(shù)(向后Euler公式)10設(shè)將方程y=f(x,y)的兩端從xn

到xn+1

求積分,得用不同的數(shù)值積分方法近似上式右端積分,可以得到計(jì)算y(xn+1)的不同的差分格式.用左矩形公式計(jì)算積分項(xiàng)再離散化，得Euler公式數(shù)值積分方法(Euler公式)11用右矩形公式計(jì)算積分項(xiàng)再離散化，得向后Euler公式數(shù)值積分方法(向后Euler公式)12用梯形法計(jì)算積分項(xiàng)再離散化，即可得如下計(jì)算公式梯形公式,單步法,隱式公式數(shù)值積分方法(梯形公式)13Taylor展開(kāi)方法(Euler公式)由此亦可得Euler公式將y(xn+1)在xn

點(diǎn)作Taylor展開(kāi),得14§1Euler方法用差分格式來(lái)近似微分方程初值問(wèn)題15例用Euler公式求解初值問(wèn)題解取步長(zhǎng)h=0.1,計(jì)算結(jié)果見(jiàn)下表計(jì)算公式為16y(xn)1.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321xn0.10.20.30.40.50.60.70.80.91.0yn1.10001.19181.27741.35821.43511.50901.58031.64981.71781.7848兩者比較可以看出Euler法的精度較差.17hhhStraightlineapproximationEuler法的幾何意義x0x1x2x3y018向后Euler法這是一個(gè)隱式公式.已知yn,上式是關(guān)于yn+1的非線(xiàn)性方程,不能從中直接求出yn+1,需要用非線(xiàn)性方程求根的方法來(lái)求解,通常用迭代法.19已知yn求yn+1的迭代法則當(dāng)hL<1時(shí),迭代函數(shù)向后Euler法計(jì)算量大,但數(shù)值穩(wěn)定性好!20Euler法的局部截?cái)嗾`差局部截?cái)嗾`差:用Euler法計(jì)算一步所產(chǎn)生的誤差,即假設(shè)第n步的計(jì)算是精確的yn=y(xn),y(xn+1)與yn+1的差,即這種誤差稱(chēng)為局部截?cái)嗾`差,簡(jiǎn)稱(chēng)截?cái)嗾`差.21當(dāng)yn=y(xn)時(shí),局部截?cái)嗾`差主項(xiàng)Euler法的局部截?cái)嗾`差22向后Euler法的局部截?cái)嗾`差定義其局部截?cái)嗾`差為向后Euler法的計(jì)算公式23向后Euler法的局部截?cái)嗾`差局部截?cái)嗾`差主項(xiàng)24在不考慮舍入誤差的情況下,稱(chēng)y(xn+1)與yn+1之差為整體截?cái)嗾`差,記為令則Euler法的整體截?cái)嗾`差25Euler法的整體截?cái)嗾`差26注意e0=0,故27所以這表明,Euler方法的整體截?cái)嗾`差與h同階.當(dāng)h0時(shí),eN

0.如果某種數(shù)值方法的局部截?cái)嗾`差為O(hp+1),則稱(chēng)它的精度是p階的,或稱(chēng)之為

p

階方法.由此可知Euler方法僅為一階方法.28梯形公式§2改進(jìn)的Euler法從直觀上看,用梯形公式計(jì)算數(shù)值積分要比矩形公式精度高.相應(yīng)地求解常微分方程的梯形公式的精度也要比Euler公式精度高.29梯形公式的局部截?cái)嗾`差局部截?cái)嗾`差主項(xiàng)故梯形公式是二階方法30則當(dāng)hL<2時(shí),梯形公式是一個(gè)隱式公式.已知yn,用迭代法計(jì)算yn+1用梯形法計(jì)算,從yn到y(tǒng)n+1要進(jìn)行迭代計(jì)算,計(jì)算量較大.實(shí)際計(jì)算常把Euler法和梯形法結(jié)合起來(lái),這就是改進(jìn)的Euler法.31先用Euler法求得一個(gè)初步的近似值，記為，稱(chēng)之為預(yù)測(cè)值，然后用它替代梯形法右端的yn+1

再直接計(jì)算fn+1,得到校正值

yn+1,這樣建立的預(yù)測(cè)－校正系統(tǒng)稱(chēng)為改進(jìn)的Euler法：預(yù)測(cè)校正改進(jìn)的Euler法32實(shí)踐表明,改進(jìn)Euler法的精度明顯高于Euler法.利用Taylor展開(kāi)可以證明,改進(jìn)Euler的局部截?cái)嗾`差為O(h3),故它是二階方法.它有下列平均化形式：33例用改進(jìn)Euler法求解解取步長(zhǎng)h=0.1,計(jì)算結(jié)果見(jiàn)下表改進(jìn)Euler法公式為34y(xn)1.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321xn0.10.20.30.40.50.60.70.80.91.0yn1.09591.18411.26621.34341.41641.48601.55251.61651.67821.7379同Euler法的計(jì)算結(jié)果比較,改進(jìn)Euler法明顯改善了精度.35改進(jìn)Euler法的幾何意義(

K1

+

K2)/2K2K1xnxn+136考察差商根據(jù)微分中值定理，存在點(diǎn)，利用所給方程y=f(x,y)

得稱(chēng)K*=f(,y())為區(qū)間[xn,xn+1]上的平均斜率,這樣只要對(duì)平均斜率K*提供一種算法，相應(yīng)地便導(dǎo)出一種計(jì)算格式.§3Runge-Kutta法Runge-Kutta法的基本思想37梯形法：改進(jìn)Euler法：其中Euler法：向后Euler法：38Runge－Kutta法的基本思想就是設(shè)法在[xn,xn+1]內(nèi)多預(yù)報(bào)幾個(gè)點(diǎn)的斜率值，然后把它們加權(quán)平均作為平均斜率，以期望構(gòu)造出更高精度的計(jì)算格式.Runge-Kutta法的基本思想39考察區(qū)間[xn,xn+1]內(nèi)一點(diǎn)xn+p=xn+ph,0<p1,用兩個(gè)點(diǎn)xn,xn+p

的斜率K1,K2的加權(quán)平均代替平均斜率K*,得如下計(jì)算格式：其中有兩個(gè)待定參數(shù),p,適當(dāng)選取它們的值，使上述格式有盡可能高的精度.二階Runge-Kutta方法40K1K2xnxn+phxn+1

=xn+h(1)

K1

+

K2加權(quán)平均斜率41若p=1/2,則該格式是二階的,故統(tǒng)稱(chēng)滿(mǎn)足這一條件的一族公式為二階Runge-Kutta方法.特別地,當(dāng)p=1,=1/2

時(shí),上述格式即為改進(jìn)的Euler法，如果取p=1/2,=1

則上述格式稱(chēng)為中點(diǎn)法.42K1K2xnxn+1/2xn+1中點(diǎn)公式

K243由此可以構(gòu)造三階Runge-Kutta方法.為了進(jìn)一步提高精度，可以考慮用三個(gè)點(diǎn)xn,xn+p,xn+q

的斜率值K1,K2,K3

加權(quán)平均得出平均斜率K*的近似值,其中三階Runge-Kutta方法44經(jīng)典的三階Runge-Kutta方法45若用四個(gè)點(diǎn)的斜率值的加權(quán)平均來(lái)近似平均斜率,可以導(dǎo)出四階Runge-Kutta方法.四階Runge-Kutta方法46經(jīng)典的四階Runge-Kutta方法47xnxn+h/2xn+hK1K2K3K4經(jīng)典的四階Runge-Kutta方法48例設(shè)取步長(zhǎng)h=0.2用四階Runge-Kutta方法求解解49y(xn)1.18321.34161.48321.61251.7321xn0.20.40.60.81.0yn1.18321.34161.48321.61241.7320同Euler法和改進(jìn)Euler法的計(jì)算結(jié)果比較,顯然以四階RK法的精度為高.四階RK法每1步計(jì)算4次函數(shù)值,改進(jìn)Euler法每1步只要計(jì)算2次函數(shù)值,但由于這里步長(zhǎng)放大了2倍,兩者所用計(jì)算量幾乎相同.取步長(zhǎng)h=0.2,計(jì)算結(jié)果見(jiàn)右表50Runge-Kutta方法的推導(dǎo)基于Taylor展開(kāi)法,因而它要求解具有較好的光滑性.如果解的光滑性差,那么,使用四階Runge-Kutta方法求得的數(shù)值解,其精度反而不如改進(jìn)的Euler法.實(shí)際計(jì)算時(shí),應(yīng)當(dāng)針對(duì)問(wèn)題的具體特點(diǎn)選擇合適的算法.51在逐步推進(jìn)的求解過(guò)程中,計(jì)算yn+1之前事實(shí)上已經(jīng)求出了一系列的近似值y0,y1,…,yn,如果充分利用前面多步的信息來(lái)預(yù)測(cè)yn+1,則可以期望會(huì)獲得較高的精度.這就是線(xiàn)性多步法的基本思想.§4線(xiàn)性多步法52線(xiàn)性二步法的一般格式局部截?cái)嗾`差53一般的線(xiàn)性多步法公式可表示為其中i,i均為常數(shù),fk=f(xk,yk)(k=n+1,n,…,nr).若1=0,則公式為顯式,若10,則公式為隱式.定義:線(xiàn)性多步法在xn+1上的局部截?cái)嗾`差定義為54基于數(shù)值積分方法非常直觀,能用它導(dǎo)出大部分的常用線(xiàn)性多步法,包括:Simpson公式

Adams顯隱公式

Milne公式構(gòu)造線(xiàn)性多步法的主要途徑是基于數(shù)值積分方法和基于Taylor展開(kāi)方法,前者可直接由方程兩端積分后利用插值求積公式得到.55Simpson公式其中將y=f(x,y)的兩端從xn1

到xn+1

求積分，得Simpson公式,兩步法,隱式公式離散得其中F(x)=f(x,y(x)),用Simpson公式計(jì)算積分項(xiàng)56Simpson公式的截?cái)嗾`差故Simpson公式是隱式二步四階方法57Adams顯式公式右端積分中用過(guò)節(jié)點(diǎn)xn,xn1,

xn2,…,

xnr的F(x)的r次插值多項(xiàng)式Lr(x)代替F(x)求積分,再離散化得(r+1)階Adams顯式公式.一階Adams顯式公式就是Euler公式58右端積分中用過(guò)節(jié)點(diǎn)xn+1,xn,

xn1,…,

xnr+1的F(x)的r次插值多項(xiàng)式Lr(x)代替F(x)求積分,再離散化得(r+1)階Adams隱式公式.Adams隱式公式一階Adams隱式公式就是向后Euler公式59二階Adams顯式公式F(x)用過(guò)節(jié)點(diǎn)

xn,

xn1的線(xiàn)性插值函數(shù)L1(x)代替60再離散化，即可得如下二階Adams顯式公式截?cái)嗾`差61二階Adams隱式公式:梯形公式再離散化，即可得如下計(jì)算公式F(x)用過(guò)節(jié)點(diǎn)

xn,

xn+1的線(xiàn)性插值函數(shù)代替62四階Adams顯式公式截?cái)嗾`差63四階Adams隱式公式截?cái)嗾`差64Milne公式F(x)用過(guò)節(jié)點(diǎn)

xn,

xn1,

xn2的2次插值函數(shù)L2(x)代替再離散化，得Milne公式65Milne公式是四階四步顯式方法,截?cái)嗾`差為66基于Taylor展開(kāi)方法能導(dǎo)出所有的常用線(xiàn)性多步法,下面以線(xiàn)性二步法為例來(lái)說(shuō)明這種方法線(xiàn)性二步法的一般格式適當(dāng)選取參數(shù)0,1,1,0,1,

使其局部截?cái)嗾`差達(dá)到最高階(使公式的階盡可能高)局部截?cái)嗾`差:假設(shè)第n步以前的計(jì)算是精確的,yi=y(xi)(in),用多步法計(jì)算一步所產(chǎn)生的誤差,即y(xn+1)與yn+1的差.67為符號(hào)簡(jiǎn)單記,把簡(jiǎn)記為由Taylor展開(kāi)式得局部截?cái)嗾`差6869要使局部截?cái)嗾`差的階高,Rn+1的前面幾項(xiàng)應(yīng)盡可能多為0,因只有5個(gè)未知數(shù),可令前5項(xiàng)為0.70其中fi=f(xi,yi),i=n

1,n,n+1.Simpson公式兩步四階隱式方法局部截?cái)嗾`差71

答案:該公式為(r+1)階Adams顯式公式適當(dāng)選取0,1,…,r

使下列公式的階盡可能高適當(dāng)選取0,1使公式的階盡可能高.例

r=1的情形局部截?cái)嗾`差72得二階Adams公式要使階盡可能高,須令局部截?cái)嗾`差73

答案:該公式為(r+1)階Adams隱式公式.適當(dāng)選取1,0,1,…,r1

使下列公式的階盡可能高74

Hamming公式四階三步隱式公式適當(dāng)選取0,2,1,0,1,

使下列公式的階盡可能高

答案:所求公式為截?cái)嗾`差75例分別用四階Adams顯式和隱式公式求解初值問(wèn)題的數(shù)值解,取h=0.1.利用線(xiàn)性多步法求解初值問(wèn)題,必須先用其他方法算出開(kāi)始幾個(gè)點(diǎn)的近似值.一般可用同階的單步法.從以上例子看到同階的Adams方法,隱式方法要比顯式方法誤差小,這可以從兩種方法的局部截?cái)嗾`差主項(xiàng)的系數(shù)大小來(lái)加以解釋.76一般地,同階的隱式法比顯式法精確,而且數(shù)值穩(wěn)定性好.但在隱式公式中,通常很難解出yn+1,需用迭代法求解,這樣又增加了計(jì)算量.因此實(shí)際計(jì)算時(shí),很少單獨(dú)使用顯式公式或隱式公式,而是將它們聯(lián)合使用:先用顯式公式求出y(xn+1)的預(yù)測(cè)值,再用隱式公式對(duì)預(yù)測(cè)值進(jìn)行校正.77仿照改進(jìn)的Euler格式的構(gòu)造方法，將顯式和隱式兩種Adams格式相匹配，可構(gòu)成下列Adams預(yù)測(cè)－校正系統(tǒng)：預(yù)測(cè)校正Adams預(yù)測(cè)-校正系統(tǒng)78預(yù)測(cè)校正Milne-Hamming預(yù)測(cè)-校正系統(tǒng)以上兩種預(yù)測(cè)-校正公式均為四階公式,其起步值通常用四階RK公式計(jì)算.79§5單步法的相容性、收斂性與穩(wěn)定性單步法的一般形式當(dāng)含有yn+1時(shí),方法是隱式的,若不含yn+1則為顯式方法.稱(chēng)(x,y,h)為增量函數(shù)顯式單步法的一般形式為80例改進(jìn)Euler法的增量函數(shù)(x,y,h)=f(x,y)例

Euler法的增量函數(shù)單步法的局部截?cái)嗾`差為p(1)階精度的單步法是指81數(shù)值解法的基本思想是,通過(guò)某種離散化手段將微分方程轉(zhuǎn)化為差分方程.轉(zhuǎn)化是否合理?即(1)當(dāng)h0時(shí),差分方程是否能無(wú)限逼近微分方程?(相容性問(wèn)題)(2)差分方程的解yn是否趨于微分方程的準(zhǔn)確解y(xn)?(收斂性問(wèn)題)82單步法可以看成是由下述方程近似而來(lái)即固定xn

,

當(dāng)h0時(shí)上式左邊極限為y(xn),如果近似是合理的,右邊應(yīng)有83定義若增量函數(shù)滿(mǎn)足(x,y,0)=f(x,y)則稱(chēng)單步法與常微分方程初值問(wèn)題是相容的.本章介紹的所有顯式單步法(RK方法)均是相容的84定義若一種數(shù)值方法對(duì)于任意固定的xn=x0+nh,當(dāng)h0(同時(shí)n

）時(shí)數(shù)值解yn趨于準(zhǔn)確解y(xn).則稱(chēng)該方法是收斂的.判斷一種數(shù)值方法的收斂性問(wèn)題本質(zhì)上就是估計(jì)該方法整體截?cái)嗾`差的問(wèn)題.例

Euler方法的整體截?cái)嗾`差為當(dāng)h0時(shí),en0,故Euler法是收斂的.收斂性85定理設(shè)單步法yn+1=yn+h(xn,yn,h)

具有p(1)階精度,其增量函數(shù)(x,y,h)關(guān)于y滿(mǎn)足Lipschitz條件又設(shè)初值是準(zhǔn)確的,即y0=y(x0),則其整體截?cái)嗾`差為y(xn)yn=O(hp)該定理意味著,滿(mǎn)足定理?xiàng)l件的單步法是收斂的只要微分方程y=f(x,y)中的函數(shù)f關(guān)于y滿(mǎn)足Lipschitz條件,則本章介紹的所有顯式單步法(RK方法)均滿(mǎn)足定理?xiàng)l件,故所有RK方法均收斂.86關(guān)于收斂性的討論有個(gè)前提，即必須假定差分方法的每一步計(jì)算都是準(zhǔn)確的(沒(méi)有舍入誤差).然而實(shí)際計(jì)算中往往由于有舍入誤差等原因而產(chǎn)生擾動(dòng)，而這些擾動(dòng)有可能“淹沒(méi)”真解,所以還要考慮穩(wěn)定性問(wèn)題.穩(wěn)定性問(wèn)題87穩(wěn)定性比較復(fù)雜.為簡(jiǎn)化討論,僅考察下列模型方程定義若一種數(shù)值方法在節(jié)點(diǎn)值yn上大小為的擾動(dòng),引起以后各節(jié)點(diǎn)值ym(m>n)上產(chǎn)生的偏差值均不超過(guò),則稱(chēng)該方法是絕對(duì)穩(wěn)定的.用某種數(shù)值方法求解上述模型方程，對(duì)固定的h,使得該方法絕對(duì)穩(wěn)定的h和的全體稱(chēng)為該方法的絕對(duì)穩(wěn)定區(qū)域.88Euler方法的絕對(duì)穩(wěn)定區(qū)域模型問(wèn)題計(jì)算格式設(shè)實(shí)際計(jì)算yn時(shí)產(chǎn)生誤差n,它使節(jié)點(diǎn)值yn+1產(chǎn)生誤差n+1,則從而為使誤差不擴(kuò)大,需89向后Euler方法的絕對(duì)穩(wěn)定區(qū)域模型問(wèn)題計(jì)算格式設(shè)實(shí)際計(jì)算yn時(shí)產(chǎn)生誤差n,它使節(jié)點(diǎn)值yn+1產(chǎn)生誤差n+1,則向后Euler方法的絕對(duì)穩(wěn)定區(qū)域90梯形公式的絕對(duì)穩(wěn)定區(qū)域模型問(wèn)題計(jì)算格式91類(lèi)似于Euler方法的討論,得梯形法的絕對(duì)穩(wěn)定區(qū)域?yàn)榧刺菪喂降慕^對(duì)穩(wěn)定區(qū)域92RK方法的絕對(duì)穩(wěn)定區(qū)域二階RK方法的絕對(duì)穩(wěn)定區(qū)域?qū)τ谀Ｐ蛦?wèn)題改進(jìn)Euler法93當(dāng)是實(shí)數(shù)時(shí),2.51h0當(dāng)是實(shí)數(shù)時(shí),2.78h0二階RK