6.2: Runge Kutta Methods (RKM) (A) 2nd Order RKM (or Improved Euler - PowerPoint PPT Presentation

6.2: Runge Kutta Methods (RKM) (A) 2nd Order RKM (or Improved Euler Method) Failure of Euler Method: Only slope on left end of interval [ t, t + h ] is used. Iteration Scheme Improvement: Given t, y ( t ), Start: y 0 , t 0 • compute slope at t s l = f ( t, y ( t )) For k = 0 to k = N : • find slope at t + h via EM t k +1 = t k + h y E = y ( t ) + hs l s l = f ( t k , y k ) s r = f ( t + h, y E ) s r = f ( t k +1 , y k + hs l ) • approximate y ( t + h ) via y k +1 = y k + h ( s l + s r ) / 2 average slope y ( t + h ) ≈ y ( t )+ h ( s l + s r ) / 2 1

Ex.: y ′ = t − y, y (0) = 0 . 5 Ex. Approximate the solution to Approximate y (1) for stepsizes y ′ = t − y, y (0) = 0 . 5 h = 1 /m, m = 1 , 2 , 4 , 8 , 16 , 32 in 0 ≤ t ≤ 1 using h = 0 . 25. Exact Value: y (1) = 0 . 551819 Start: y 0 = 0 . 5, t 0 = 0 Error: E ( h ) = | y (1) − y m | t 1 = 0 . 25 E ( h ) h y m s l = t 0 − y 0 = − 0 . 5 1 0 . 75 0 . 198181 s r = t 1 − ( y 0 + hs l ) = − 0 . 125 1 / 2 0 . 585938 0 . 034118 y 1 = y 0 + h ( s l + s r ) / 2 = 0 . 4219 1 / 4 0 . 558794 0 . 006974 1 / 8 0 . 553400 0 . 001581 t 2 = 0 . 5 1 / 16 0 . 552196 0 . 000377 = t 1 − y 1 = − 0 . 1719 s l 1 / 32 0 . 551911 0 . 000092 s r = t 2 − ( y 1 + hs l ) = 0 . 1211 y 2 = y 1 + h ( s l + s r ) / 2 = 0 . 4155 E ( h/ 2) ≈ E ( h ) / 4 ⇒ E ( h ) ≈ C h 2 t 3 = 0 . 75 Theorem: There ∃ C > 0 s.t. s l = t 2 − y 2 = 0 . 0845 = t 3 − ( y 2 + hs l ) = 0 . 3134 s r E ( h ) ≤ C h 2 y 3 = y 2 + h ( s l + s r ) / 2 = 0 . 4653 t 4 = 1 (2nd order RKM is second s l = t 3 − y 3 = 0 . 0845 order method) s r = t 4 − ( y 3 + hs l ) = 0 . 3134 = y 3 + h ( s l + s r ) / 2 = 0 . 4653 y 4 2

(B) 4th Order RKM Ex.: y ′ = t − y, y (0) = 0 . 5, y (1) ≈ y m Idea: Given t and y = y ( t ), compute slopes s 1 , s 2 , s 3 , s 4 at m = 1 , 2 , 4 , 8 , 16 , 32 , h = 1 /m four carefully chosen points s.t. Exact Value: y (1) = 0 . 551819162 error is minimized. Error: E ( h ) = | y (1) − y m | h y m E ( h ) Approximation: 1 0 . 5625 0 . 010680838 y ( t + h ) ≈ y + h 6( s 1 +2 s 2 +2 s 3 + s 4 ) 1 / 2 0 . 552256266 0 . 000437105 1 / 4 0 . 551841299 0 . 000022137 1 / 8 0 . 551820408 0 . 000001246 Iteration k → k + 1 : 1 / 16 0 . 551819236 0 . 000000074 1 / 32 0 . 551819166 0 . 000000005 s 1 = f ( t k , y k ) s 2 = f ( t k + h/ 2 , y k + hs 1 / 2) E ( h/ 2) ≈ E ( h ) / 16 ⇒ E ( h ) ≈ C h 4 s 3 = f ( t k + h/ 2 , y k + hs 2 / 2) Theorem: There ∃ C > 0 s.t. s 4 = f ( t k + h, y k + hs 3 ) E ( h ) ≤ C h 4 y k +1 = y k + (4th order RKM is fourth h 6( s 1 + 2 s 2 + 2 s 3 + s 4 ) order method) t k +1 = t k + h 3

Error Comparison Least square fit of E ( h ) EM: E(h) ≈ 0.3412 h EM Ex. y ′ = t − y, y (0) = 0 . 5 RK2: E(h) ≈ 0.1348 h 2 0.15 y (1) ≈ y m → E ( h ) = | y (1) − y m | 0.1 h = 1 /m, m = 1 , 2 , 4 , 8 , 16 , 32 E(h) h EM RKM2 RKM4 0.05 RK2 1 0 . 5518 0 . 198181 0 . 010680838 1 / 2 0 . 1768 0 . 034118 0 . 000437105 0 1 / 4 0 . 0772 0 . 006974 0 . 000022137 0 0.1 0.2 0.3 0.4 0.5 h 1 / 8 0 . 0364 0 . 001581 0 . 000001246 − 4 x 10 1 / 16 0 . 0177 0 . 000377 0 . 000000074 RK4: E(h) ≈ 0.0070 h 4 1 / 32 0 . 0087 0 . 000092 0 . 000000005 4 3 E ( h ) for E(h) EM: Euler Method 2 RKM2: 2nd order RKM 1 RKM4: 4th order RKM 0 0 0.1 0.2 0.3 0.4 0.5 h 4

Worked out Examples from Exercises Ex. 3: y ′ = ty , y (0) = 1. Compute five RK2-iterates for h = 0 . 1. Arrange computation and results in a table. k t k y k s l s r h h ( s l + s r ) / 2 0 0 1 0 0 . 1 0 . 1 0 . 005 1 0 . 1 1 . 0050 0 . 1005 0 . 2030 0 . 1 0 . 0152 2 0 . 2 1 . 0202 0 . 2040 0 . 3122 0 . 1 0 . 0258 3 0 . 3 1 . 0460 0 . 3138 0 . 4309 0 . 1 0 . 0372 4 0 . 4 1 . 0832 0 . 4333 0 . 5633 0 . 1 0 . 0498 5 0 . 5 1 . 1331 0 . 5665 0 . 7138 0 . 1 0 . 0640 5

Ex. 7: z ′ + z = cos x , z (0) = 1 (i) Compute RK2-approximations in 0 ≤ x ≤ 1 for h = 0 . 2, h = 0 . 1, h = 0 . 05. (ii) Find exact solution (iii) Plot exact solution as curve and RK2 approximations as points. (i) In Matlab, the RK2 approximation for h = 0 . 2 is computed and stored � x in arrays x0 2 , z0 2 via e − x + e ξ cos( ξ ) d ξ z ( x ) = h=0.2; 0 (cos x + sin x + e − x ) / 2 m=1/h;x=0;z=1; = xv=x;zv=z; for k=1:m (iii) Plot: sl=cos(x)-z; (see CN Sec. 6.1 for commands) sr=cos(x+h)-(z+sl*h); z=z+h*(sl+sr)/2;zv=[zv z]; 1 x=x+h;xv=[xv x]; 0.98 end 0.96 o: h=0.2 x0_2=xv;z0_2=zv; 0.94 z *: h=0.1 0.92 Analogously for h = 0 . 1 and h = 0 . 05 +: h=0.05 0.9 (arrays x0 1 , z0 1 and x0 05 , z0 05 ). 0.88 0.86 (ii) Variation of Parameter: 0 0.2 0.4 0.6 0.8 1 x h = − z ⇒ z h ( x ) = e − x z ′ 6

Ex. 7a: z ′ + z = cos x , z (0) = 1 (i) Compute RK4-approximation in 0 ≤ x ≤ 1 for h = 0 . 2. (iii) Plot exact solution as curve and RK4 approximation as points. (iii) Matlab plot commands: (i) RK4 approximation for h = 0 . 2 is x=linspace(0,1,100); computed and stored in arrays xv , zv : z=1/2*(cos(x)+sin(x)+exp(-x)); h=0.2; plot(xv,zv,’ko’,x,z,’k’), m=1/h;x=0;z=1; xlabel(’x’),ylabel(’z’), xv=x;zv=z; axis([0 1 0.86 1]) for k=1:m s1=cos(x)-z; 1 s2=cos(x+h/2)-(z+s1*h/2); 0.98 s3=cos(x+h/2)-(z+s2*h/2); 0.96 s4=cos(x+h)-(z+s3*h); 0.94 z=z+h*(s1+2*s2+2*s3+s4)/6; z 0.92 zv=[zv z]; 0.9 x=x+h;xv=[xv x]; 0.88 end 0.86 0 0.2 0.4 0.6 0.8 1 x 7