MIT OpenCourseWare http://ocw.mit.edu 16.346 Astrodynamics Fall 2008 For information about citing these materials or our Terms of Use, visit: http://ocw.mit.edu/terms.Lecture 33 Runge-Kutta Method of Numerical Integration #12 Numerical Integration dr dx dxi= v = y = yidt = dt dt i =1, 2, 3 dv ⇒dy ⇐⇒dyi = g(r) = f(x) = fi(xi )dt dt dt x(t0)= x0 y(t0)= y0 f(x0)= f0 Taylor Series Expansion dx h2 d2x h2x(t + h)= x0 + h ��+ � ���� �+ ··· = x fdt 2 0=tdt0 + hy0 + + �t0 2! t=t0 2! ···dy y(t + h)= y0 + h��h2 d2y ���h2 � + � + ··· = y + hf f dt 0 0 + + t=t0 2! dt2 t=t0 2! 0···h2 ∂f dx = y0 + hf0 + x ���� + t=t0··· 2! ∂ dt= F0y0 First Order Method Second Order Method � �� � x = x2 2 30 + hy1 0 + O(h ) x = x0 + hy0 + h f0 + O(h2) y = y h20 + hf0 + O( ) y = y0 + hf0 + 1 h22F0y + 30 O(h ) How to Avoid Calculating the Matrix F0 ∂f f(x0 + δx)= f0 + ���� δx + O[(δx)2]∂x t=t0 or for any constant p ∂f f(x0 + hpy0)= f��1�+ O(h20 + hpy0 = f0 + hpF0y0)∂x t=t0 or [f(x� 0 + hpy0) − f0]= hF0y0 + O(h2) p Hence, the second-order method is equivalent to 12 3� � x = x0 + hy0 + 2h f(x0)+ O(h ) 1 1 h f(x0 + hpy0)+ O(h3)+ h 1 − f(x0)+y = y0 2p 2p Choose p = 1 and note that f(x0) − f(x0 + hpy0)= O(h). Therefore, we have 2 x = x0 + hy0 + 12h2f(x0 + 12hy0)+ O(h3) y = y0 + hf(x0 + 12hy0)+ O(h3) with only one evaluation of the function f(x) for x = x0 + 12hy0. 16.346 Astrodynamics Lecture 330 0 1 0 0 21Equate corresponding coefficients of α0 and α� � � � 1 a = 1 1 1 b a p 2 = 1 = ⇒ = = 1 b 2 =1 p 2 In summary: x = x1 0 + hy0 + h22k + O(h3) y = y0 + hk + O(h3) k = f(x0 + 1 h2y0) Taylor Expansion using Indicial Notation and Summation Convention dxi h2 d2xi h3 d3xi h4 d4xi xi(t + h)=x i+ h + + + +dt 2! dt2 3! dt3 4! dt4 ··· 2 3 i 4 2 i = i + i h hdf hd fx hy + fi + + + 2! 3! dt 4! dt2 ···3 32 3 i j 4 i j= i + i h+ i h� ∂f dxh� d ∂fx hy fdt� dx+ + + 2! 3! ∂xj dt 4! ∂xjdt j=1 j=1 � ···h2 h3 h4 = x i + fi+ j y + hyi + fij(fiyjy k+ fi jjk jf )+2! 3! 4! ··· Taylor Expansion of a Vector Function of a Vector fi(xi + δi)= fi+ fiδj+ 1 fiδjδk+ 1 fiδjjkδkδj2 6 jk+ ··· 16.346 Astrodynamics Lecture 33 ���� ���� � � Formal Derivation of the Second-Order Method Choose a, b, p so that x = x0 + hy0 + h2 ak + O(h3) x = x0 + hy0 + 12h2α0 + O(h3) y = y0 + hbk + O(h3) ≡ y = y0 + h(α0 + 1hα1)+ O(h3)2k = f(x0 + hpy0) k = α0 + hpα1 + O(h2) where α0 = f0 and α1 = f0= F0y0 . Expand f(x0 + hpy0) in a Taylor series: ∂f ∂f δx + O[(δx)2] = f(x0 + hpy0)= f0 + +O(h2)f(x + δx)= f0 + hpy0⇒ ∂x ∂x t=t0 ��t=t0 = α0 + hpα1Then x = x0 + hy0 + h2 a(α0 + hpα1)+ O(h3) x = x0 + hy0 + 12h2α0 + O(h3) y = y + hb(α + hpα )+ O(h3) ≡ y = y + h(α + 1hα )+ O(h3)Nystr¨om’s Third Order Method Evert Johannes Nystr¨om 1925 x = x0 + hy0 + h2(a0k0 + a1k )+ O(h41) y = y40 + h (b0k0 + b1k1)+ O(h ) k0 = f(x0 + hp0y0) k1 = f(x0 + hp2 1y0 + h q1k0) Series Expansion Using Indicial Notation xi (t + h)= x i + hyi + 1h2fi + 13 h fi yj2 6j+ O(h4) i i i 12 i j 13 i j k i j 4� �� � h f + h (f + fy(t + h)= y + hf + 2jy6jky yjf )+ O(h ) k0 i = fi(x i + hp0y i)= fi + fihp0yj + 1fi (hp0yj)(hp0y k)+ O(h3)� �� � j2jkδi 0 k1 i = fi(x i + hp1y i + h2 q1k0 i )= fi[x i + hp1y i + h2 q1fi + O(h3)] δi = fi + fi(hp1yj + h2 fj)+ 1fi (hp1yj)(hp1y k)+ O(h3)1 jq12jkSeries Expansion Using Vector Notation x = x0 + hy0 + h2(12α0 + 16hα1)+ O(h4) y = y0 + h[α0 + 12hα1 + 16h2(α2 + β2)] + O(h4) k0 = f(x0 + hp0y0)= α0 + hp0α1 + 12h2 p02α2 + O(h3) k1 = f(x0 + hp1y0 + h2 q1k0)= α0 + hp1α1 + h2(21 p12α2 + q1β2)+ O(h3) Determine a0 , a1 , b0 , b1 , p0 , p1 , q1 h2(a0k0 + a1k1) ≡ h2(1α0 + 1hα1) to terms of order h3 2 6h(b0k0 + b1k1) ≡ h[α0 + 21hα1 + 61h2(α2 + β2)] to terms of order h3 Condition Equations � �� � � 1 �⎡ 1 1 ⎤� b0 � ⎡ 1 ⎤ (α ) 1 1 a0 = 21 ⎣ p0 p1 ⎦ b1 = ⎣ 12 ⎦ p0 p1 a16 2 21 p0 p1 3 (β ) q1b1 = 12 ( 13 ) 16.346 Astrodynamics Lecture 330 1 1 1 p b + p b = 1 (1 p0)b0 +(1 p1)b1 = a2 0 + a1 = 2 0 0 1 1 2 = − −⇒ p0(1 − p0)b b⇐⇒= 1 0 + p1(1 p1)1 = 1 p0a0 + p1a1 p 2 0b0 + p 2 1b1 1 = 6 6 3 −Hence a0 =(1 − p0)b0 and a1 =(1 − p1)b1 . Next, for consistency, b0 + b1 =1 ���1 1 1 p0b0 + p1b1 1 = 2 3 = �p 2= 1 �0b0 + p 2 ⇐⇒ D1b1 3��p p1 p��Expand the determinan�0 1 �� 2 20 p 2113 �=0 t � D3 = ��1 V��p0 p���1 1 = � �����andermonde determinant 1 (p1 11 2 1 2 21 �� − p0� )[ − (p0+ p1)+p0p1] = V2 L3(p0,p1)Complete � p 0 p 13 solution � �3 2Constrain��t function � �=0 �� � 1 p1 b2 1 p2 0 1 p1 p0 0 = −b1 = −q1 =−a0 =(1 − p0)b0 a1 =(1 − p1)bp1 0 − p1 p1 − p0 6 1 − p20 with p0 and p1 chosen arbitrarily subject to L (p ,p )= 1− 1 3 0 1(p3 20 + p1)+ p0p1 =0 Nystr¨om’s Third Order Algorithm with p0 =0 x = x0 + hy0 + 1h24(k0 + k1)+ O(h4) y = y0 + 1h O4(k0 +3k1)+ (h4) k0 = f(x0) k = …