3 MATLAB Runge-Kutta Butcher 3. Taylor Taylor y(x 0 + h) = y(x 0 ) + h y (x 0 ) + h! y (x 0 ) + Taylor 3. Euler, Runge-Kutta Adams Implicit Euler, Implicit Runge-Kutta Gear y n+ y n (n+ ) y n+ y n+ y n+ 3.3 d dt d dt d dt [A] = k[a] [B] = k[a] k[b] [C] = k[b] 4
d dt A B C = 0 0 0 0 0 A B C d dt x = Ax, A 0 B 0 C 0 = mm 0 mm 0 mm AP = PΛ P AP = Λ P A = ΛP P d dt x = P Ax d dt (P x) = Λ(P x) d dt P x = 0 0 0 0 0 0 0 P x P x = = C C exp( t) C 3 exp( t) 0 0 0 A B C A + B + C = C A + B = C exp( t) A = C 3 exp( t) 5
A 0 B 0 C 0 = C 3 = C = C = mm 0 mm 0 mm A + B + C = A + B = exp( t) A = exp( t) A B C = exp( t) exp( t) exp( t) exp( t) + exp( t) mm 3.4 3.4. : Euler Euler Taylor y(x 0 + h) = y(x 0 ) + h y (x 0 ) + h! y (x 0 ) + Euler y(x 0 + h) = y(x 0 ) + h y (x 0 ) Euler Taylor t S S(t) t t + t S S(t + t) Euler S(t + t) = S(t) + ts (t) S (t) S ds dt 6
S(t + t) = S(t) + t ds dt ds dt ds = k S dt S(t + t) = S(t) t k S S(t) = S(t + t) S(t) = t k S Euler t S t k S 3.4.. d[glucose] = k[glucose] dt t =.0sec t = 0 5.0 sec k = 0.sec t = 0 sec [Glucose]=.0mM. () t = 5.0 sec ε = n X analytical i Xi numerical n X analytical i i= 3. t.0 sec 4. 7
3.4.3 : Euler 3.4.4 3.5 α 0 y n + α y n + + α k y n k = h(β 0 f n + β f n + + β m f n m ) y n t n y y n k = y(t n kh) f = y 3.5. Euler RK 3.5. : Adams-Bashforth Gear Adams y n = α y n + h(β f n + β f n ) α, β Taylor y(t n ) = α y(t n ) + h(β f(t n, y(t n )) + β f(t n, y(t n ))) t n = t h Taylor y(t n ) α { y(t n ) h d dt y(t n) + h + hβ { f(t n ) h d dt f(t n) + h + hβ { f(t n ) h d dt f(t n) + 4h } d dt y(t n) d dt f(t n) } d dt f(t n) } 8
h h 0 : y(t n ) = α y(t n ) h : 0 = α hy (t n ) + hβ y (t n ) + hβ y (t n ) h : 0 = h α y (t n ) h β y (t n ) h β y (t n ) α, β = α 0 = α + β + β 0 = α β β Adams-Bashforth α =, β = 3, β = Adams-Bashforth ds dt S = t ( 3 y n = y n + h f n + ) f n = k S ( 3 k S t=t t ) k S t=t t 3.5.3 : Adams-Moulton 3.5.4 (predictor-corrector) Adams-Bashforth Adams-Moulton MATLAB ode3 ( ) ( ) Backward Differential Formula(ode5s) 9
3.6 3.6. Runge-Kutta 3.7 Runge-Kutta Runge-Kutta 3.7. LM LM RK Adams 3.7. s Runge-Kutta k = f(x 0, y 0 ) k = f(x 0 + c h, y 0 + ha k ) k 3 = f(x 0 + c 3 h, y 0 + h(a 3 k + a 3 k )). k s = f(x 0 + c s h, y 0 + h(a s k + a s,s k s )) y = y 0 + h(b k + + b s k s ) 0
c = a c 3 = a 3 + a 3. c s = a s + + a s,s s Runge-Kutta 3.7.3 Butcher s Runge-Kutta Butcher k 0 k c a k 3 c 3 a 3 a 3...... k s c s a s a s a s,s y b b b s b s 3.7.4 Runge-Kutta (RK) Runge-Kutta k = f(x 0, y 0 ) k = f(x 0 + h, y 0 + hk ) y = y 0 + hk 3.7.5 Runge-Kutta y(x 0 + h) y Kh p+ Runge-Kutta p Runge-Kutta Taylor p
3.7.6 Taylor Euler y = y 0 + hf ( x 0 + h, y 0 + h ) f 0 Taylor f ( ( ) x 0 + h, y 0 + hf 0) = h f(x0, y 0 ) + + hf x 0 f(x y 0, y 0 ) ( + h + hf! x 0 y) f(x0, y 0 ) + = f(x 0, y 0 ) + h (f x + ff y ) (x 0, y 0 ) + ( h 4 + h f x 0 x + h f y 4 0 y ) f(x 0, y 0 ) + = f(x 0, y 0 ) + h (f x + ff y ) (x 0, y 0 ) + h (f 8 xx + ff xy + f f yy ) (x 0, y 0 ) + y = y 0 + hf(x 0, y 0 ) + h (f x + f y f) (x 0, y 0 ) + h3 (f 8 xx + f xy f + f yy f ) (x 0, y 0 ) + Taylor y(x 0 + h) = y 0 + hy (x 0, y 0 ) + h y (x 0, y 0 ) + h3 6 y (x 0, y 0 ) + = y 0 + hf(x 0, y 0 ) + h x + f y f)(x 0, y 0 ) + h3 6 (f xx + f xy f + f yy f + f x f y + f y f)(x 0, y 0 ) + y(x 0 + h) y = h3 4 {(f xx + f xy f + f yy f + 4(f x f y + f y f)}(x 0, y 0 ) + y(x 0 + h) y Kh 3 Runge-Kutta
3.7.7 dy dx = f(x, y) df = f x df dx = f x + f y = y y = f x + f y f dx + f y dy dy dx dy = y x dy dx dx + y y dy = y x + y y = dy dx x (f x + f y f) + y (f x + f y f) dy dx = {f xx + (f xy f + f y f x )} + {f xy + (f yy f + fy )}f y = f xx + f xy f + f yy f + f x f y + f y f 3.8 3.8. p(ˆp) p ˆp ( : Dormand-Prince 5(4), Bogacki-Shampine 3() ode45, ode3 ) Dormand-Prince Runge-Kutta(5,4) 3.8. RK3() RK3 3 3 3
k = f(x n, y n ) () k = f(x n + c h, y n + ha k ) () k 3 = f(x n + c h, y n + h(a 3 k + a 3 k )) (3) y n+ = y n + h(b k + b k + b 3 k 3 ) (4) k 3 h 3 k (4) h h ( ) k = f(x n, y n ) + h c! x k f y ( + h c! x k y) f + O(h 3 ) = f(x n, y n ) + h(c f x + a ff y ) + h (c f xx + c a ff xy + a f f yy ) + O(h 3 ) k 3 k 3 = f(x n, y n ) + h! + h! ( c 3 ) f + (a x 3k + a 3 k ) y ( c 3 + (a x 3k + a 3 k ) y) f + O(h 3 ) = f(x n, y n ) + h(c 3 f x + a 3 ff y + a 3 k f y ) + h {c 3f xx + c 3 (a 3 f + a 3 k )f xy + (a 3 f + a 3 k ) f yy + O(h 3 ) k f a 3 k a 3 k = a 3 { f + h(c f x + a ff y ) + O(h ) } a 3 k = a 3 (f + O(h)) h k 3 k 3 = f(x n, y n ) + h[c 3 f x + a 3 ff y + a 3 {f + h(c f x + a ff y ) + O(h )} f y ] + h [c 3f xx + c 3 {a 3 f + a 3 (f + O(h))}f xy + {a 3 f + a 3 (f + O(h))} f yy ] RK h RK 4
y(x + h) = y(x) + h dy dx + h d y! dx + h3 d 3 y 3! dx + 3 O(h4 ) y(x + h) y(x) = hf + h (f x + f y f) + h3 6 (f xx + f xy f + f yy f + f y f x + fy f) + O(h 4 ) d y dx = f x + f y f df dx = f x + f df y dx RK h RK : hf = h(b + b + b 3 )f : 3 : h (f x + f y f) = b h (c f x + a ff y ) + b 3 h (c 3 f x + a 3 ff y + a 3 ff y ) h 3 6 (f xx + f xy f + f yy f + f y f x + fy f) = b h 3 (c f xx + c a ff xy + a f f yy ) + b 3 h 3 a 3 (c f x + a ff y )f y + b 3h 3 {c 3f xx h +c 3 a 3 ff xy + c 3 a 3 ff xy +(a 3 + a 3 ) f f yy } b + b + b 3 = h h f x : f x = b h c f x + b 3 h c 3 f x = b c + b 3 c 3 h ff y : f yf = b h a ff y + b 3 h (a 3 + a 3 )ff y = b a + b 3 (a 3 + a 3 ) h 3 5
h 3 f xx : 6 f xx = b h 3 c f xx + b 3h 3 c 3f xx = 6 b c + b 3c 3 h 3 ff xy : 6 (f xyf) = b h 3 (c a ff xy ) + b 3h 3 (c 3 a 3 ff xy + c 3 a 3 ff xy ) = b c a + b 3 c 3 (a 3 + a 3 ) 3 f yy f h 3 : 6 f yyf = b h 3 a f f yy + b 3h 3 (a 3 + a 3 ) f f yy = 6 b a + b 3(a 3 + a 3 ) h 3 f y f x : 6 f yf x = b 3 h 3 a 3 c f x f y = b 3 a 3 c 6 fy h 3 f : 6 f y f = b 3 h 3 a 3 a ffy = b 3 a 3 a 6 c = a c 3 = a 3 + a 3 b + b + b 3 = b c + b 3 c 3 = b c + b 3 c 3 = 3 b 3 a 3 c = 6 8 (a, a 3, a 3, b, b, b 3, c, c 3 ) 6 c = u, c 3 = v { b u + b 3 v = b u + b 3 v = 3 { b u + b 3 uv = u b u + b 3 v = 3 6
(uv v )b 3 = u 3 b 3 = = u 3 uv v 3u 6v(u v) { b uv + b 3 v = v b u + b 3 v = 3 (uv u )b = v 3 3v b = 6u(v u) b 3 a 3 c = 6 3u 6v(u v) a 3u = 6 a 3 = v(u v) u(3u ) b = b b 3 a = c = u a 3 = c 3 a 3 = v a 3 (RK) Runge-Kutta b, b, b 3, ŷ n+ = y n + h(ˆb k + + ˆb s k s ) Runge- Kutta 7
ε tol h ( ) RK3() ˆb 3 = 0 y i ŷ i ε tol b u + b 3 v = ˆb 3 = 0 ˆb = ˆb ˆb u = ˆb = = u c ˆb = c (RK) 0 0 c =, a =, ˆb = 0, ˆb = RK3 0 c 3 a 3 a 3 b b b 3 RK3 Ralston (REF.FIXME) 0 3 4 0 3 4 9 3 Ralston a 3 = 0, a 3 = 3, b 4 =, b 9 =, b 3 3 = Ralston RK3 4 9, c 3 = 3 4 (REF.FIXME) 4 9 8
3.8.3 y i ŷ i ε tol err = n n ( yi ŷ i ε i= tol ) err err Ch q+ Ch q+ opt err ( h h opt ) q+ h opt h( q+ err h opt. err h. err h opt h 3.9 3.9. Gear (BDF: Backward Differential Formula) α 0 y n + α y n + + α k y n k = hf n Adams Taylor Gear Numerical initial value problems pp.4 stiffly stable method 9
L.H.S. α 0 y(t n ) + α { y(t n ) h d dt y(t n) + h + α { y(t n ) h d dt y(t n) + 4h R.H.S. = h d dt y(t n) } d dt y(t n) d dt y(t n) } h h 0 : α 0 y(t n ) + α y(t n ) + α y(t n ) = 0 h : α hy (t n ) hα y (t n ) = hy (t n ) h : α h y (t n ) + α h y (t n ) = 0 α 0 + α + α = 0 α α = + α = 0 Gear α 0 = 3, α =, α = Gear ds dt 3 y n y n + y n = hf n = k S 3 S t S t t + S t t = t k S t 3.9. (stiff equation) (stiffness ratio) = > 0 4 (stiff) 30
A BDF A (stiffly stable) 6 BDF 3.0 Further reading Gear, Springer s yellow book http://www.scholarpedia.org/article/backward_differentiation_formulas by Gear 3