計(jì)算方法4常微分方程數(shù)值解法

1、常微分方程的數(shù)值解法Numerical Solutions to Ordinary Differential Equations 對象對象一階常微分方程初值問題一階常微分方程初值問題: 00y)x(y)y,x( fdxdy一階常微分方程組初值問題一階常微分方程組初值問題: 000n0n202101n21nnn2122n2111y)x(y,y)x(y,y)x(y)y,y,y,x(fy)y,y,y,x(fy)y,y,y,x(fy高階常微分方程初值問題高階常微分方程初值問題: )n(00)n(0000)n()n(y)x(y,y)x(y,y)x(y)y,y,y,x( fy(4.1)一階常微分方程初值問

2、題一階常微分方程初值問題:實(shí)際工程技術(shù)、生產(chǎn)、科研上會出現(xiàn)大量的微分方程問實(shí)際工程技術(shù)、生產(chǎn)、科研上會出現(xiàn)大量的微分方程問題很難得到其解析解，有的甚至無法用解析表達(dá)式來表題很難得到其解析解，有的甚至無法用解析表達(dá)式來表示，因此只能依賴于數(shù)值方法去獲得微分方程的數(shù)值解。示，因此只能依賴于數(shù)值方法去獲得微分方程的數(shù)值解。 00y)x(y)y,x( fdxdy常數(shù)常數(shù)稱為稱為有唯一解有唯一解則初值問題則初值問題使得使得條件：條件：上連續(xù)，且滿足上連續(xù)，且滿足在在若若定理定理LipschitzL,)1.4(|yy|L|)y,x( f)y,x( f|,y,y,LLipschitzb,a)y,x( f21

3、2121 ab x0 x1 x2 . xn-1 xn用數(shù)值方法，求得用數(shù)值方法，求得y(x)在每個節(jié)點(diǎn)在每個節(jié)點(diǎn)xk 的值的值y(xk ) 的近似的近似值，用值，用yk 表示，即表示，即yk y(xk)，這樣，這樣y0, y1,.,yn稱為微分稱為微分方程的數(shù)值解方程的數(shù)值解求求y(x)求求y0, y1,.,yn 微分方程的微分方程的數(shù)值解法數(shù)值解法：不求不求y=y(x)的精確表達(dá)式的精確表達(dá)式,而求離散點(diǎn)而求離散點(diǎn)x0,x1,xn處的處的函數(shù)值函數(shù)值設(shè)設(shè)(4.1) 的解的解y(x)的存在區(qū)間是的存在區(qū)間是a,b，初始點(diǎn)，初始點(diǎn)x0=a，取取a,b內(nèi)的一系列節(jié)點(diǎn)內(nèi)的一系列節(jié)點(diǎn)x0, x1,.,

4、xn。a= x0 x1 xn =b，一般采用等距步長，一般采用等距步長思路思路計(jì)算過程計(jì)算過程：方法方法：采用步進(jìn)式和遞推法：采用步進(jìn)式和遞推法將將a,bn等分等分, a= x0 x14)階龍格庫塔法。階龍格庫塔法。 m4時：計(jì)算量太大，精確度不一定提高，有時會降低。時：計(jì)算量太大，精確度不一定提高，有時會降低。 Gill公式公式 )hK)221(hK22y,hx( fK)hK)221(hK212y,2hx( fK)K2hy,2hx( fK)y,x( fK)KK)22(K)22(K(6hyy32nn421nn31nn2nn14321n1n節(jié)省存儲單元節(jié)省存儲單元控制舍入誤差控制舍入誤差 對于經(jīng)

5、典的四階對于經(jīng)典的四階Runge-Kutta法給出如下算法：法給出如下算法：算法算法4.2求解：求解： dy/dx=f(x,y) axb y (a)=y0 Step 1: 輸入輸入a,b,y0 及及NStep 2: (b-a)/N=h,a=x,y0=yStep 3: 輸出輸出 (x,y)Step 4: For I=1 T0 Nhf(x,y)=k1hf(x+h/2,y+ k1/2)= k2hf(x+h/2,y+k2/2)=k3hf(x+h,yk3)=k4y+(k1+2k2+2k3+k4)/6=yx+h=x輸出輸出(x,y)Step 5 : 停止停止例例4.3用四階經(jīng)典用四階經(jīng)典RungeKutt

6、a方法解初值問題：方法解初值問題： 102yyxydxdy2 . 0h(1)求求 ， 1y2 . 0, 1, 000hyx12,000001yxyyxfK9181818. 022222,2100101002KhyhxKhyKhyhxfK9086375. 022222,2200202003KhyhxKhyKhyhxfK8432399. 02,300303004hKyhxhKyhKyhxfK1832293. 126432101KKKKhyy 1832160. 14 . 11211xxy51103 . 1e(2)求求 ， 2y10.2,0.2xh1832293. 11y8451714. 02,111

7、111yxyyxfK211122,0.7946656hhKf xyK311222,0.7874945hhKf xyK4113,0.7440375Kf xh yhK21123421.34168036hyyKKKK3416408. 11222xxy52100 . 4e自適應(yīng)龍格庫塔自適應(yīng)龍格庫塔法法用戶提出問題I ：問題：問題：如何判斷：如何判斷|y(xn)-yn| 精度值精度值y(xn)未知。未知。 ：如何取如何取h=?解解：如用：如用p階龍格庫塔法計(jì)算，局部截斷誤差為階龍格庫塔法計(jì)算，局部截斷誤差為O(hp+1) y=f(x,y) y(x0)=y0 要求誤差要求誤差10-8求數(shù)值解求數(shù)值解 x

8、n h/2 xn+1h如如 xn-xn+1 令令 yn=y(xn) yn+1(h) 則則 y(xn+1)-yn+1(h)chp+1步長折半步長折半xnxn+h/2xn+1分兩步計(jì)算分兩步計(jì)算y(xn+1)的近似值的近似值yn+1(h/2)。則則 y(xn+1)-yn+1(h/2)2c(h/2)p+1 h2n+1( )n+1(h)pn+1n+1y(x)-y1y(x)-y2h2h2nn+1+1( )(h)n+1pn+112y -yyy1(x)-定理：對于問題I若用P階龍格庫塔法計(jì)算y(xn+1)在步長折半前后的近似值分別為yn+1(h), yn+1(h/2)則有誤差公式 ph2n+1h2n+112

9、( )n( )(h)n+111+y(x|)-yyy |-=|注：注：10 誤差的事后估計(jì)法誤差的事后估計(jì)法 20 停機(jī)準(zhǔn)則：停機(jī)準(zhǔn)則： (可保證可保證|y(xn+1)-yn+1(h/2)|） 將步長折半反復(fù)計(jì)算，直至將步長折半反復(fù)計(jì)算，直至為止，為止， 取取h為最后一次的步長為最后一次的步長， yn+1為最后一次計(jì)算的結(jié)果為最后一次計(jì)算的結(jié)果。 Else if (為止，最后為止，最后一次運(yùn)算的前一次計(jì)算結(jié)果即為所需。一次運(yùn)算的前一次計(jì)算結(jié)果即為所需。 4.5 收斂與穩(wěn)定性收斂與穩(wěn)定性 對于常微分方程初值問題對于常微分方程初值問題(4.1)求求y=y(x)-求求y(xn)- 求求yn xn=x0

10、+nh離散化離散化某種公式某種公式例：例：Euler法，改進(jìn)的法，改進(jìn)的Euler法，龍格庫塔法都是單步法。法，龍格庫塔法都是單步法。數(shù)值解法：數(shù)值解法： 00y = fx, yyx= y單步法：計(jì)算單步法：計(jì)算yn+1時只用到前一步的結(jié)果時只用到前一步的結(jié)果yn。顯式單步法：顯式單步法： yn+1=yn+h(xn,yn,h) (x,y,h)為增量函數(shù)，它依賴于f,僅是xn,yn, h的函數(shù)Def:若某數(shù)值方法對于任意固定的若某數(shù)值方法對于任意固定的xn=x0+nh,當(dāng)當(dāng) h-0(n-)時，時， yn-y(xn)，則稱該法是則稱該法是收斂收斂的。的。0(0)(limnhnny xy即即 (xn

11、=x0+nh為固定值為固定值)改進(jìn)改進(jìn)Euler法的收斂性法的收斂性 判定改進(jìn)判定改進(jìn)Euler法的收斂性：法的收斂性： yn+1=yn+hf(xn,yn)+f(xn+h,yn+hf(xn,yn)/2 則則(x,y,h)=f(x,y)+f(x+h,y+hf(x,y)/2若：若：yo=y(x0) _|f(x,y)-f(x,y)| L|y-y|f(x,y)關(guān)于關(guān)于y滿足李滿足李-條件條件:12_ |(x,y,h)- (x,y,h)|f(x,y)-f(x,y)|+|f(x+h,)-f(x+y+hf(x,y)y+hh,f(x,y)|則：則： 12_L|y- y|+L|y+hf(x,y)- y-hf(x

12、,y)|_21L|y- y|+L|y- y|+L hhL(1+L)2|y- y|=|y- y|2限定限定hh0，則則_0_| (x,y,h)- (x,y,hhL(1+L)2)|y-y|即即 (x,y,h) 滿足李普希茲條件滿足李普希茲條件,由由Th改進(jìn)的改進(jìn)的Euler法收斂法收斂同樣可驗(yàn)證，若同樣可驗(yàn)證，若f(x,y)關(guān)于關(guān)于y滿足李普希茲條件滿足李普希茲條件,且且 y0=y(x0)時，時， Euler法，標(biāo)準(zhǔn)四階龍格庫法也收斂。法，標(biāo)準(zhǔn)四階龍格庫法也收斂。4.5.2 穩(wěn)定性穩(wěn)定性在討論收斂性時，一般認(rèn)可：數(shù)值方法本身計(jì)算在討論收斂性時，一般認(rèn)可：數(shù)值方法本身計(jì)算過程是準(zhǔn)確的。實(shí)際上，并非如

13、此：過程是準(zhǔn)確的。實(shí)際上，并非如此： 初始值初始值y0有誤差有誤差0y0-y(x0). 計(jì)算的每一步有舍入誤差。計(jì)算的每一步有舍入誤差。這些初始誤差在計(jì)算過程的傳播中，是逐步衰減這些初始誤差在計(jì)算過程的傳播中，是逐步衰減的，還是惡性增長，這就是穩(wěn)定性問題。的，還是惡性增長，這就是穩(wěn)定性問題。Def4.1若一種數(shù)值方法在節(jié)點(diǎn)若一種數(shù)值方法在節(jié)點(diǎn)xn處的數(shù)值解處的數(shù)值解yn的擾動的擾動n0，而在以后的各節(jié)點(diǎn)值而在以后的各節(jié)點(diǎn)值ym(mn)上有擾動上有擾動m.當(dāng)當(dāng)|m|n|，(m=n+1,n+2,)，則稱該數(shù)值算法則稱該數(shù)值算法是是穩(wěn)定穩(wěn)定的。的。分析某算法的數(shù)值穩(wěn)定性，通?？疾炷Ｐ头匠谭治瞿乘惴ǖ?/p>

14、數(shù)值穩(wěn)定性，通?？疾炷Ｐ头匠?y=y ， (0)Def4.1 ： 設(shè)在節(jié)點(diǎn)設(shè)在節(jié)點(diǎn)xn處用數(shù)值法得到的理想數(shù)值解為處用數(shù)值法得到的理想數(shù)值解為yn，而實(shí)際計(jì)算得到的近似值為，而實(shí)際計(jì)算得到的近似值為 ，稱，稱 為第為第n步數(shù)值解的擾動步數(shù)值解的擾動nynnnyy Euler法的穩(wěn)定性法的穩(wěn)定性Euler法：yn+1=yn+hf(xn,yn)考察模型方程考察模型方程 y=y，(0) 即即yn+1=(1+h)yn假設(shè)在節(jié)點(diǎn)值假設(shè)在節(jié)點(diǎn)值yn上有擾動上有擾動n，在，在yn+1上有擾動上有擾動n+1，且且n+1僅由僅由n引起引起(計(jì)算過程不再引進(jìn)新的誤差計(jì)算過程不再引進(jìn)新的誤差)Euler法穩(wěn)定法穩(wěn)定

15、Euler法的穩(wěn)定的條件為法的穩(wěn)定的條件為 0h-2/nn1n1n1ny)h1(y)h1(yy nnn1n)h1()yy)(h1( 1h1 n1n 隱式隱式Euler法穩(wěn)定性法穩(wěn)定性隱式隱式Euler法法：yn+1=yn+hf(xn+1,yn+1)對于模型方程對于模型方程 y=y( yn+1=yn/(1-h)yn+1的擾動的擾動0 上式均成立上式均成立，隱式隱式Euler法無條件穩(wěn)定法無條件穩(wěn)定h1h1yh1yyynnn1n1n1n 1h11Euler 法法穩(wěn)穩(wěn)定定隱隱式式梯形公式穩(wěn)定性梯形公式穩(wěn)定性梯形公式梯形公式 yn+1=yn+hf(xn,yn)+f(xn+1,yn+1)/2，設(shè)模型方程

16、為設(shè)模型方程為y=y(0)則則 yn+1=yn+hyn+yn+1/2 n+1nh1 +2y=yh1 -20時恒成立時恒成立 梯形公式恒穩(wěn)定。梯形公式恒穩(wěn)定。 yn+1的擾動的擾動n1n2h12h1 12h12h1 梯梯形形公公式式穩(wěn)穩(wěn)定定 00y = fx, yyx= yTaylor公式公式公式公式單步否單步否顯式否顯式否精度精度 n+1nnny=y +hf(x ,y )單步單步顯式顯式一階一階n+1nn+1n+1y=y +hf(x,y)單步單步隱式隱式一階一階n+1n-1nny=y+2hf(x ,y )二步二步顯式顯式二階二階數(shù)值積分法數(shù)值積分法hn+1nnnn+1n+12y=y + f(x

17、 ,y )+f(x,y)單步單步隱式隱式二階二階11n+1n12221nn2nn1y= y +k +kk = hf(x ,y )k = hf(x +h,y +k )單步單步顯式顯式二階二階RungeKutta方法方法局部截斷誤差局部截斷誤差整體截斷誤差整體截斷誤差收斂與穩(wěn)定性收斂與穩(wěn)定性4.6 多步法多步法某一步解和前若干步解的值有關(guān)某一步解和前若干步解的值有關(guān)1N,m,1mn)y,x( fbhyaymj1nmj1n1m0jm0jjmj1nj1n m步線性多步法的一般形式步線性多步法的一般形式其中其中a0,am-1, b0,bm為和為和h無關(guān)的常數(shù)無關(guān)的常數(shù)初值為初值為y0, y1, ,ym-

18、1若若bm0，方法是隱式的方法是隱式的(implicit)否則是顯式的否則是顯式的(explicit)(4.7)AdamsBashforth四階方法四階方法1N,4,3n)y,x( f9)y,x( f37)y,x( f59)y,x( f55(24hyy,y,y,y,y3n3n2n2n1n1nnnn1n3322110 AdamsMoulton四階方法四階方法1N,3,2n)y,x( f)y,x( f5)y,x( f19)y,x( f9(24hyy,y,y,y2n2n1n1nnn1n1nn1n22110 多步法的截斷誤差多步法的截斷誤差1N,m,1mn)x(y,x( fb)x(ya)x(y(h1)

19、h(7.4m0jmj1nmj1nj1m0jmj1nj1n1n ）式式的的截截斷斷誤誤差差為為定定義義（4.7 常微分方程組和高階微分方程的數(shù)值解法常微分方程組和高階微分方程的數(shù)值解法形如形如 0000z)x(z),z,y,x(gzy)x(y),z,y,x( fy歐拉方法的計(jì)算公式歐拉方法的計(jì)算公式 00nnnn1n00nnnn1nz)x(z),z,y,x(hgzzy)x(y),z,y,x(hfyy隱式歐拉方法的計(jì)算公式隱式歐拉方法的計(jì)算公式 001n1n1nn1n001n1n1nn1nz)x(z),z,y,x(hgzzy)x(y),z,y,x(hfyy梯形公式梯形公式 )z,y,x(g)z,y

20、,x(g(2hzz)z,y,x( f)z,y,x( f (2hyy1n1n1nnnnn1n1n1n1nnnnn1nAdams公式公式R-K方法方法對于高階微分方程，總可以化成方程組的形式，對于高階微分方程，總可以化成方程組的形式，例如例如 0000y)x(y,y)x(y)y,y,x( fy可以化成可以化成 0000z)x(z,y)x(y)z,y,x( fzzy例例6.0)0(y,4.0)0(y,1x0,xsiney2y2yx2 令令y(x)=z，則有則有 z2y2xsinezzyx2用用4階階RK方法，有方法，有6.0z,4.0y,0 x000 04.0)z2y2xsine(h)z,y,x(h

21、fk06.0hz)z,y,x(hfk000 x200022,1000011,10 70324764475.0)2kz,2ky,2hx(hfk062.0)2kz(h)2kz,2ky,2hx(hfk2,101,10022,22,102,101,10011,2 70315240923.0)2kz(2)2ky(2)05.0 xsin(e(hk80616283223.0)2kz(hk2,201,200)05.0 x(22,32,201,30 80217863729.0)kz(2)ky(2)1.0 xsin(e(hk40631524092.0)kz(hk2,301,300)1.0 x(22,42,301,

22、40 6/)kkkk(zz6/)kkkk(yy2,42,32,22,1011,41,31,21,101 4.8 常微分方程邊值問題的數(shù)值解法常微分方程邊值問題的數(shù)值解法形如形如bxa),y,y,x( fy 第一類邊界條件第一類邊界條件 )b(y,)a(y(4.8)第二類邊界條件第二類邊界條件 )b(y,)a(y第三類邊界條件第三類邊界條件 )b(y,)a(y定理定理 假定問題（假定問題（4.8）式中的函數(shù)）式中的函數(shù)f在集合在集合 y,y,bxa)y,y,x(D上連續(xù)，且偏導(dǎo)數(shù)上連續(xù)，且偏導(dǎo)數(shù)fy和和fy也在也在D上連續(xù)。若有上連續(xù)。若有M)y,y,x(fM) ii(0)y,y,x(f ,D)

23、y,y,x() i (yy 使使得得存存在在常常數(shù)數(shù)對對于于所所有有則邊值問題有唯一解。則邊值問題有唯一解。例例 邊值問題邊值問題0)2(y)1(y,2x1,0ysineyxy 有有ysine)y,y,x( fxy 因此邊值問題有唯一解。因此邊值問題有唯一解。1ycos)y,y,x(f,0 xe)y,y,x(fyxyy 定理定理 若線性邊值問題若線性邊值問題 )b(y,)a(y,bxa),x( ry)x(qy)x(py滿足滿足0)x(qb,a) ii()x( r)x(q),x(pb,a) i ( 上上在在連連續(xù)續(xù)和和上上在在則邊值問題有唯一解。則邊值問題有唯一解。打靶法打靶法（Shooting Method）)x(y)b(y)b(y)x(y)x(y2211 其中，其中，y1(x)是初值問題是初值問題的解的解0)a(y,)a(y,bxa),x( ry)x(qy)x(py 的解，的解，y2(x)是初值問題是初值問題1)a(y,0)a(y,bxa,y)x(qy)x(py 有限差分法有限差分法(Finite-Differe