數(shù)值分析慶課件群共享第六章版

1、Introduction to Numerical AnalysisChapter 6: Numerical Methods for ODEs慶本人所有2012 年 5 月 14 日慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日1 / 83常微分微分描述自然界的連續(xù)變化。常微分描述的是未知函數(shù)同其導數(shù)之間的。例如(牛頓運動定律)：F = ma F (t, y (t ), dy (t )/dt ) = md 2y (t )/dt 2.一般形式：F (x, y (x ), y 0 (x ), y 00 (x ), &

2、#183; · · , y (n)(x ) = 0.慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日2 / 83常微分初值問題初值問題（IVP）：(y 0 = f (x, y ),y (a) = .a x b,f (x, y ), f (x, y )假設連續(xù)。y在此假設下，上述IVP的等價形式：唯一解y (x ) 且在 a, b 上充分光滑。Zxfy (x ) = +(t, y (t )dt.a慶 (本人所有)Introduction to Numerical Analysis Chapte20

3、12 年 5 月 14 日3 / 83為什么需要數(shù)值解？很多，如分離變量法，因子法，級數(shù)解法等。但遺憾的是，實際得到的大部分微分都不能使用這些理論上的，通常只能得到解的某種近似。即使得到解，也可能因形式復雜而不具有實際操作意義。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日4 / 83數(shù)值解將 a, b n 等分，記h = (b a)/n,稱 h 為步長。xi = a + ih(i = 0, 1, · · · , n),所謂數(shù)值解，是求初值問題的解 y (x ) 在離散點 xi 處的

4、近似值 yi。計算計算法。yi +1yi +1時，如果只用到的值 yi ，稱為單步法。時需用到前 r 步的值 yi , yi1, · · · , yir+1， r 步方慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日5 / 83數(shù)值解和解的區(qū)別雖然大多數(shù)情況下，兩種可能滿足要求，但是：解可以給出任意時刻的解的計算公式；數(shù)值解只能給出解函數(shù)在離散點上的近似值列表；無窮數(shù)值解為有限上的函數(shù)；上的離散。慶 (本人所有)Introduction to Numerical Analysis Ch

5、apte2012 年 5 月 14 日6 / 83數(shù)值解的過程離散化：a x0 < x1 < x2 < · · · < xn b。構(gòu)造近似列表：數(shù)值求解 ODE 的通常是構(gòu)造遞推公式。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日7 / 83在數(shù)值求解之前分析原析解，的信息：是否滿足解的解的形式是否足夠好等等。唯一性條件，是否有解分析解的穩(wěn)定性。這是因為在由一個離散點的值計算下一個點的值時，通常都會產(chǎn)生一定的誤差。解的穩(wěn)定性決定了這類誤差隨著時間的增加放大或者

6、縮小。分析數(shù)值算法的時，還必須關(guān)注算法的穩(wěn)定性和收斂性。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日8 / 83Euler 公式（1）問題重述：(y 0 = f (x, y ),y (a) = .a x b,將兩邊在 xi , xi+1ZZxi +1y 0xi +1f(x )dx(x, y (x )dx,=xixi即Zxi +1fy (xy (x(x, y (x )dx.i +1) =i ) +xi慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5

7、 月 14 日9 / 83Euler 公式（2）應用左矩形公式近似右端得(1)y (x) = y (x ) + hf (x , y (x ) + R,i +1iiii +1其中1 df (x, y (x )h2 = 1 y 00(i )h2,R(1)i (xi , xi=1), +i +12dx2x =iR(1)上式中忽略有i +1y (xi+1) y (xi ) + hf (xi , y (xi ) yi + hf (xi , yi ),慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日10 / 83公式（3）Eu

8、ler得 Euler公式：yi+1 = yi + hf (xi , yi ),i = 0, 1, · · · , n 1.此公式稱為向前 Euler 公式。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日11 / 83ffff圖示：Introduction to Numerical Analysis Chapte2012 年 5 月 14 日12 / 83ff慶 (本人所有)Eulerfffor y 0 = ymethody. . . . . . . . . . . . . .y0 =

9、 y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .y0 . . . . . . .tttttt01234Euler 公式的幾何解釋Euler 公式的其它導出方式向前 Euler公式導出方式還有很多，如： Taylor級數(shù)法、有限差分法、多項式插值法等。Taylor級數(shù)法：有限差分法：慶 (本人所有)Introduct

10、ion to Numerical Analysis Chapte2012 年 5 月 14 日13 / 83單步法的一般形式Euler 公式是一個單步顯式公式。一般的單步顯式公式為：(yi+1 = yi + h(xi , yi , h),y0 = .(x, y, h) 稱為增量函數(shù)。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日14 / 83ffEuler公式的另一個例子：Euler method for y 0 = yy.y0. . . . . .0y = y. . . . . . . . . . . . .

11、. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .t.ttttt01234Introduction to Numerical Analysis Chapte2012 年 5 月 14 日15 / 83ffffff慶 (本人所有)Euler 公式的另一個例子ODE數(shù)值解中的誤差舍入誤差。截斷誤差。在絕大多數(shù)情形下，截斷誤差占據(jù)著主導地位。以后略舍入誤差。忽但是要注意，當步長 h 太小時，舍入誤差將占據(jù)主要位置。截斷誤差

12、又分為局部截斷誤差和整體誤差兩種。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日16 / 83局部誤差（1）稱Ri+1 = y (xi+1) y (xi ) + h(xi , y (xi ), h)為單步顯式公式在xi +1點的局部截斷誤差。由上述定義，向前 Euler 公式的局部截斷誤差為：Ri+1 = y (xi+1) y (xi ) + hf (xi , y (xi ),1=h2y 00(i ),i (xi , xi 1).+2慶 (本人所有)Introduction to Numerical Analys

13、is Chapte2012 年 5 月 14 日17 / 83局部誤差（2）局部截斷誤差指的是單步產(chǎn)生的誤差。它不考慮誤差的。在假設求解精確的情況下，計算算法產(chǎn)生的誤差。不同問題的局部截斷誤差對整體截斷誤差的影響也會不同。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日18 / 83ff圖示：error for y 0 = yy. . . . . . .global error. . . . . . . . . . . . .0. . . .y = y. . . . . . . . . . . . . . . .

14、. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .y.0 . . . . .tt0t1t2t3t4Introduction to Numerical Analysis Chapte2012 年 5 月 14 日19 / 83慶 (本人所有)局部誤差（3）ff圖示：error for y 0 = yy.y.0 . . . . . . .y0 = y. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15、. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .global errort. .t0t1t2t3t4Introduction to Numerical Analysis Chapte2012 年 5 月 14 日20 / 83慶 (本人所有)ff局部誤差（4）后退（Backward）Euler公式（1）用右矩形公式近似得(2)y (x) = y (x ) + hf (x, y (x) + R,i +1ii +1i +1i +1其中2h2 y 00h df (x, y (x ) R(2)= 2(i ),i (x

16、i , xi +1 = 2i +1). dxx =i從而有y (xi+1) y (xi ) + hf (xi+1, y (xi+1) yi + hf (xi+1, yi+1),慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日21 / 83后退（Backward）Euler 公式（2）后退 Euler 公式為：yi+1 = yi + hf (xi+1, yi+1),i = 0, 1, · · · , n 1.單步隱式公式。一般的單步隱式公式為yi+1 = yi + h(xi , yi ,

17、 yi+1, h),y0 = .其中 (xi , yi , yi+1, h) 稱為增量函數(shù)。i = 0, 1, · · · , n 1慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日22 / 83后退（Backward）Euler 公式的截斷誤差由定義，后退 Euler 公式的局部截斷誤差為：Ri +1= y (xi+1) y (xi ) hf (xi+1, y (xi+1)h2= 2 y (i ),i (xi , xi+1).00慶 (本人所有)Introduction to Num

18、erical Analysis Chapte2012 年 5 月 14 日23 / 83梯形格式（1）用梯形公式近似得y (xi) = y (x ) + h f (x , y (x ) + f (x, y (x) + R(3) ,+1iiii +1i +1i +12其中32h d f (x, y (x )1 y 000 R(3)h3= 12(i ),i (xi , x= 12 i +1).i +1dx 2x =i略去R(3) 得i +1hy (xi 1) y (xi ) + 2 f (xi , y (xi ) + f (xi 1, y (xi 1)+h yi + 2 f (xi , yi )

19、+ f (xi+1, yi+1),慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日24 / 83梯形格式（2）梯形公式y(tǒng)i 1 = yi + h f (xi , yi ) + f (xi 1, yi 1),i = 0, 1. · · · , n 1.+2它是一個單步隱式公式。由單步隱式公式局部截斷誤差的定義得梯形公式的局部截斷誤差為：= y (xi 1) y (xi ) + 2 f (xi , y (xi ) + f (xi 1, y (xi 1) hRi+1+ 1 12y 0003(

20、 )h ,i (xi , xi 1).= i+慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日25 / 83梯形格式的缺點梯形格式y(tǒng)i 1 = yi + h f (xi , yi ) + f (xi1, yi 1),i = 0, 1. · · · , n 1.+2是一個隱式算法。隱式算法導致每一步都要反解一個。該還很有可能是非線性的。這導致工作量加大，而且會引入新的誤差。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5

21、月 14 日26 / 83校值：y (p)= yi + hf (xi , yi ).i +1代替梯形格式中右邊的yi +1。校正：= y + h f (x , y ) + f (x, y (p) ).yi+1iiii +1i +12慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日27 / 83改進的Euler 公式即是：y (p)= yi + hf (xi , yi )i +1y (c)(p)i +1= y + hf (x, y)ii +1i +1(p)(c)1y= (y+ y)i +12i +1i +1稱為改進的

22、Euler 公式。單步顯式公式。也可以表示為：= y + h f (x , y ) + f (xyi, y + hf (x , y ).+1iiii +1iii2慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日28 / 83改進的 Euler 公式的截斷誤差（1）其局部截斷誤差為：Ri +1= y (xi+1)y (x ) +f (x , y (x ) + f (x 1, y (xi ) + hf (xi , y (xi )i . hhiiii +2可用兩種求上面的局部截斷誤差。慶 (本人所有)Introducti

23、on to Numerical Analysis Chapte2012 年 5 月 14 日29 / 83改進的Euler 公式的截斷誤差（2）法一：Ri +1h= y (xi+1) y (xi ) 2 f (xi , y (xi ) + f (xi+1, y (xi+1)+ 2 f (xi+1, y (xi+1) f (xi+1, y (xi ) + hf (xi , y (xi )h 1 12 1 12h f 2(x, )i+1iy+1y 0003( )h +y (x 1) y (xi ) hf (xi , y (xi )= = ii +h 1 f 2 2(x, )3i+1iy+1y 00

24、2y 000( )h +( )hii 1 121 f 4(x, )i+1iy+13( ) +( ) h ,i , (xi , xi 1).y 000y 00= iii+慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日30 / 83改進的Euler 公式的截斷誤差（3）法二：將 y (xi+1) 在 xi 進行處 Taylor 展開：得到： 1 12y 003y 0004y (xi) = y (x ) + hy 0(x ) +h(x ) +h(x ) + O(h ),+1iiii2!3!慶 (本人所有)Introdu

25、ction to Numerical Analysis Chapte2012 年 5 月 14 日31 / 83改進的 Euler 公式的截斷誤差（4）將 f (xi+1, y (xi ) + hf (xi , y (xi ) 在 (xi , y (xi )開：得到：進行處Taylor展f (xi+1, y (xi ) + hf (xi , y (xi ) f f= f (xi , y (xi ) + hx + hf (xi , y (xi ) y 12 2f 2f2 fy 22223h x 2 + 2h f (xi , y (xi ) xy + h (f (xi , y (xi )+ O(h

26、 )+ 2! 1fy2 y 0003y 0hy 00y 00=(x ) +(x ) + h(x ) (x )+ O(h ).iiii2慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日32 / 83改進的 Euler 公式的截斷誤差（5）將上面 2 式代入誤差式并利用 y (x ) 及其導數(shù)和f (x, y (x )的，得到： 1 12y 003y 0004Ri= y (x ) + hy 0(x ) +h(x ) +h(x ) + O(h )+1iiii2!3!1y (x ) hy 0(x )ii2 fy112 y

27、0003 h(x ) +(x ) + hy 0hy 00(x ) y 00(x )+ O(h)iiii22 f 1 12134y 000y 00(x ) +(x )h + O(h ).= iiy4其截斷誤差是 O(h3) 的。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日33 / 83ODE 數(shù)值解公式的階如果一個求解公式的局部截斷誤差為 O(hp+1)，則稱該公式是p 階的，或具有 p 階精度。根據(jù)這定義，Euler 公式、后退的 Euler 公式是 1 階的，梯形公式和改進的 Euler 公式是 2 階的。慶

28、 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日34 / 83回顧：改進的 Euler 公式= y + h f (x , y ) + f (x梯形格式：yi, y)。+1iiii +1i +12值：y (p)= yi + hf (xi , yi )。i +1= y + h f (x , y ) + f (x, y (p) )。校正算法：yi幾何解釋：+1iiii +1i +12慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日35 / 83

29、平均斜率（1）由微分中值定理，y (xi+1) y (xi ) = y 0(xi + h), (0, 1),h利用得：y (xi+1) = y (xi ) + hf (xi + h, y (xi + h),稱 f (xi + h, y (xi + h) 為 y (x ) 在 xi , xi+1 上的平均斜率。記為 k 。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日36 / 83平均斜率（2）記k1 = f (xi , yi ),k2 = f (xi+1, yi + hk1),若用若用改進k1 近似 k ，則得

30、1 階Euler公式；k1 + k2近似 k ，則得 2 階改進的Euler公式。2Euler 公式精度高的：多取了一個點的斜率值，對平均斜率估計較好。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日37 / 83算法的想法（1）R-KZxi +1f已知：y (xyi +(x, y (x )dx。i +1) =xi只要算的精確些，yi+1 自然會精確嗎？問題：如何提高增加節(jié)點。格式的精度呢？慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14

31、日38 / 83算法的想法（2）R-K假設Zrxi +1f j(t, y (t )dt hc f (xi + h, y (xi + h).jjxij =1如何處理 y (xi + j h) 呢？一個思路是：利用已知估計未知！慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日39 / 83顯式 R-K 算法（1）一般的 r級 Runge-Kutta為：r= y + h kyi +1ij jj =1k1 = f (xi , yi )k = f x + h, y + h j 1µ k,l =1j = 2, 3,

32、· · · , r .jijijl l慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日40 / 83顯式 R-K 算法（2）選擇參數(shù) j , j , µjl 使其具有一定的階數(shù)。局部截斷誤差rRi+1 = y (xi+1) y (xi ) h j Kj ,j =1其中K1 = f (xi , y (xi ),j = f x!j 1i + h, y (xi ) + hµ K,Kj = 2 3, , · · · , r ,jjl ll =1

33、慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日41 / 83顯式 R-K 算法（3）將截斷誤差展開為 h 的冪級數(shù)Ri+1 = c0 + c1h + · · · + cphp + cp+1hp+1 + · · ·選擇參數(shù) j , j , µjl ，使得 c0 = c1 = · · · = cp = 0，而cp+1 6= 0，則公式是 p 階的。慶 (本人所有)Introduction to Numerical An

34、alysis Chapte2012 年 5 月 14 日42 / 83的階R-KRi +1 = O(hp+1)時，得到r 級 p 階R-K算法。r = 1 時，取 b1 = 1 得到是 Euler 算法。如何確定算法的階呢？工具：Taylor 展開。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日43 / 832 階 Runge-Kutta 公式（1）2階Runge-Kutta 公式一般形式為yi+1 = yi + h(1k1 + 2k2)k1 = f (xi , yi )k2 = f (xi + 2h, yi +

35、 hµ21k1)其局部截斷誤差是：Ri+1 = y (xi+1) y (xi ) h(1K1 + 2K2)K1 = f (xi , y (xi )K2 = f (xi + 2h, y (xi ) + hµ21K1)慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日44 / 832 階 Runge-Kutta 公式（2）利用 Taylor 展開等，得到：K1 = y 0(xi ),K2 = f (xi , y (xi ) + 2h f + hµ21y 0(xi ) fxy 2 2f 2f2

36、1 fy 22 0023+ 2 (2h) x 2 + 22µ21h y (x )+ (µ hy (x )+ O(hi21iy21 12y 003y 0004y (xi) = y (x ) + hy 0(x ) + h(x ) +h(x ) + O(h )+1iiii23! 1ff124= y (x ) + hy 0(x ) + h+(x )+(x ) + O(h )y 0y 000iii2xy6慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日45 / 832 階 Runge-Kutta 公式（3

37、）將上面3 式 代入局部截斷誤差得，Ri +1= h(1 1 2)y 0(xi ) f1 fy122+h + µy 0(x )2 22 21ix2" 2 2f 2fxy113000+h6 y (xi ) 2 2 (2)+ 2 µ y (x )02 21ix2# f2024+(µ21y (xi )+ O(h ).y 2慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日46 / 832 階 Runge-Kutta 公式（4）要使R-k算法具有 2 階精度，則1 1 2 = 0121

38、2 = 02 2 µ= 02 21即1 = 1 22 = 1 22µ21 = 1 .22其中2 6= 0。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日47 / 832 階 Runge-Kutta 公式（5）從而得到一類2階 Runge-Kutta 公式y(tǒng)i+1 = yi + h(1 2)k1 + 2k2k1 = f (xi , yi ) 1 1 k = fx +h, y +hk.2ii12222慶 (本人所有)Introduction to Numerical Analysis Chapte

39、2012 年 5 月 14 日48 / 832 階 Runge-Kutta 公式（6）2 = 1 ，得改進的 Euler 公式。當當22 = 1，得變形的 Euler 公式：yi+1 = yi + hk2k1 = f (xi , yi ) 11k = fx + h, y + hk.2ii1222 = 3 ，則得若4= y + h k1 + 3k2)yi+1i4k1 = f (xi , yi ) 22k = fx + h, y + hk.2ii133慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日49 / 83其他

40、R-K 算法及隱式 R-K 算法類似的利用上述構(gòu)造Runge-Kutta 公式?？梢缘玫?3 階或4階等高階有隱式的 R-k 算法，如：(yi+1 = yi + hk1k1 = f (xi + 1 h, yi + 1 k1h)22隱式算法具有更好的數(shù)值穩(wěn)定性。在求解 stiff組時經(jīng)常要使用隱式的 R-K算法。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日50 / 83整體誤差hnn設 y (x )是微分的解，y是用某種數(shù)值得ii =1ii =1到的近似解。則稱hE (h) =max |y (x ) y|ii1i n為該的整體截斷誤差。如果lim E (h) = 0,h0則稱該收斂。慶 (本人所有)Introduction to Numerical Analysis Chapte2012 年 5 月 14 日51 / 83單步法的收斂性（1）考慮單步顯式公式(yi+1 = yi + h(xi , yi , h),y0 = .i = 0, 1, ·