University of Hyogo 8 8 1 d x(t) =f(t, x(t)), dt (1) x(t 0 ) =x 0 () t n = t 0 + n t x x n n x n x 0 x i i = 0,..., n 1 x n x(t) 1 1.1 1 1 1 0 θ 1 θ x n x n 1 t = θf(t n 1, x n 1 ) + (1 θ)f(t n, x n ) () Euler Euler θ = 1 Euler Euler Euler x n x n 1 t = f(t n 1, x n 1 ) (4) Euler Euler θ = 1 Euler x n x n 1 t = f(t n, x n ) (5) x n x n 1
Trapezoidal Rule θ = 1/ x n x n 1 t = 1 [f(t n 1, x n 1 ) + f(t n, x n )] (6) 1. dx dt = f(t, x), x(t 0) = x 0 (7) t [t 0, t 1 ] f t [t 0, t 1 ], x x 0 b M t 1 t 0 < b/m t, x, x Lipshitz f(t, x) f(t, x ) L x x (8) N t = t 1 t 0 Euler (x(t)) 1 C x(t n ) x n C t (9) Proof. Euler x n+1 = x n + tf(t n, x n ). (10) (7) t n t n+1 tn+1 x(t n+1 ) = x(t n ) + f(t, x)dt = x(t n ) + tf(t n, x(t n )) + O( t ), (11) t n (f(t) F (t) F (t + t) F (t) = tf(t) + O( t )) x(t n+1 ) x n+1 = x(t n ) x n + t[f(t n, x(t n )) f(t n, x n )] + O( t ). (1) Lipshitz x(t n+1 ) x n+1 x(t n ) x n + t f(t n, x(t n ) f(t n, x n ) + A t. (1) f(t n, x(t n ) f(t n, x n ) L x(t n ) x n (14) e n x(t n ) x n e n+1 e n + tl e n + A t = (1 + tl) e n + A t (15) e 0 = 0 n 1 e n A t (1 + tl) i = A t L [(1 + tl)n 1] (16)
(1 + tl) N = exp(l(t 1 t 0 )) (17) x(t n ) x n A t L [exp(l(t 1 t 0 )) 1] (18) 1. (local truncation error) t n 1 x n 1 = x(t n 1 ) t n t n x(t n ) x n e = x(t n ) x n (19) Scheme L(x n 1,, x n k )x n = 0 (0) e = Lx(t n ) (1) Scheme e = O( t m ) or e A t m () Scheme m m ( ) m 1 Euler x(t n+1 ) Taylor x(t n+1 ) = x(t n + t) = x(t n ) + t dx dt + t d x dt + O( t ) () dx/dt = f e Euler = O( t ). (4) Euler 1 Euler 1 Fractional Step Method 1 t n 1 x n 1 f(t, x(t)) x n
.1 Runge-Kutta R- (stage) Runge-Kutta R k (r) = f(t n 1 + ta r, x n 1 + t b rs k (s) ) (r = 1,,, R) (5) s=1 R x n = x n 1 + t c r k (r). (6) r=1 a r, b rs, c r R a r = b rs, r = 1,,, R (7) ( X = s=1 ) ( t, F (X) = x(t) ) 1 (8) f(t, x(t)) (5),(6) (consistency condition) (5),(6) R c i = 1. (9) i=1 1. b rs (s r). b rs (s > r) k (1), k (),. R R 1 4 5 6 7 8 9 10 R 10 m 1 4 4 5 6 6 7 7 m R R Euler a 1 = b 11 = 0, c 1 = 1 Euler a 1 = b 11 = 1, c 1 = 1 1 Runge-Kutta 1. a = b 1 = 1/, c = 1 k (1) = f(t n 1, x n 1 ) (0) k () = f(t n 1 + t, x n 1 + t k(1) ) (1) x n = x n 1 + tk () () 4
. a = b 1 = 1, c 1 = c = 1/ 1 k () = f(t n 1 + t, x n 1 + t f(t n 1, x n 1 )) = f + t = f + t x(t n ) Taylor k (1) = f(t n 1, x n 1 ) () k () = f(t n 1 + t, x n 1 + tk (1) ) (4) x n = x n 1 + t (k(1) + k () ) (5) f t + t f f x + 1 ( t f 4 t + t 4 f f x t + t 4 f f x ( f t + f f ) ( + t f x 8 t + f f x t + f f x ) + O( t ) ) + O( t ) (6) x(t n ) = x(t n 1 ) + t dx dt + t d x dt + t d x 6 dt + O( t4 ) (7) dx = f dt d x dt = df dt = f t + f f x d x dx = d f dt = d ( f dt t + f f ) x = ( f t t + f f ) ( f + x t + f f x = f f + (f + ) f t t x + f f t x + f ) f x + f ( f x t + f f ) x ( f x (8) (9) (40) (41) ) + f f x (4) e rk = x(t n+1 ) x n+1 = O( t ) (4) k (1) + k () = f + t 1 O( t ) ( f t + f f ) ( + t ) f x 4 t + f f x t + f f x + O( t ) (44) 4 4 Runge-Kutta Runge-Kutta 5
(a 1 = 0, a = b 1 = 1/, a = b = 1/, a 4 = b 4 = 1, c 1 = c = c = c 4 = 1/) 4 k (1) = f(t n 1, x n 1 ) (45) k () = f(t n 1 + t, x n 1 + t k(1) ) (46) k () = f(t n 1 + t, x n 1 + t k() ) (47) k (4) = f(t n 1 + t, x n 1 + tk () ) (48) x n = x n 1 + t 6 [k(1) + k () + k () + k (4) ] (49) e rk4 = x(t n ) x n = t5 ( ) (50) 880 ( )4 4 Runge-Kutta Gill (Runge-Kutta-GIll ) do i = 1, 4 u = f(x) x = x + t (P i u + Q i v) v = R i u + S i v end do ( = ) u, v Runge-Kutta ((45)-(48) ) Runge-Kutta P 1 = b 1 P = b, P 1 + Q R 1 = b 1 + Q R 1 = b 1 (5) P = b 4, P + Q R = b + Q R = b 4, (5) (51) P 1 + Q R 1 + Q S R 1 = b 1 + Q S R 1 = b 41 (54) P 4 = c 4, P + Q 4 R = b 4 + Q 4 R = c, (55) P + Q R + Q 4 S R = b 4 + Q 4 S R = c, (56) P 1 + Q R 1 + Q S R 1 + Q 4 S S R 1 = b 41 + Q 4 S S R 1 = c 1 (57) Runge-Kutta (51)-(57) (b 41 b 1 )(c b 4 ) = (b 4 b )(c 1 b 41 ). (58) 6
a, b, c a 1 = 0, a = 1, b 1 = 1 a = 1, b 1 = c 1 6c, b = 1 6c a 4 = 1, b 41 = 0, b 4 = 1 c, b 4 = c (59) c 1 = 1 6, c = c, c, c 4 = 1 6 c 4 Runge-Kutta a 4 = 1 c = + 6 a 1 = 0, a = 1, b 1 = 1 a = 1, b 1 = 1+, b = a 4 = 1, b 41 = 0, b 4 =, b 4 = + c 1 = 1 6, c = 6, c = + 6, c 4 = 1 6 (60) a, b, c 8 (Q, Q, Q 4, R 1, R, R, S, S ) (51)-(57) 5 R 1 = 1 Q = b 1 b 1 Q = b 4, Q 4 = c 4 P 1 = 1 P = P = + P 4 = 1 6 Q 1 = 0 Q = Q = + Q 4 = 1 R 1 = 1 R = R = + R 4 = 0 S 1 = 0 S = 4 S = +4 S 4 = 0 (61) i = 1 v Q 1 = S 1 = 0 v R 4 = S 4 = 0 Q 1, S 1 Linear Multistep Method ( x f(t, x(t)) x n k β 0 = 0 (8) (6) α i x n i = t β i f(t n i, x n i ) (α 0 0) (6) α i t n i = t β i (6) i i α i x n i = t β i f(t n i, x n i ) (64) i i 7
t n t n 1 = t t n α i = t (β i + iα i ) (65) i i t n, t α i = 0 (66) iα i = (66),(67) (consistency condition) k β i (67) γ i x n i = ϕ n, n = n 0, n 0 + 1, (68) γ n γ 0 0, γ k 0 x n0, x n0+1, {x n } γ i x n i = 0, n = n 0, n 0 + 1, (69) {ˆx n } (68) {ψ n } x n = ˆx n + ψ n ϕ n k γ i 0 ψ n = ϕ/ γ i (70) (69) {x n,t }, t = 1,, K a 1 x n,1 + a x n, + + a K x n,k = 0, n = n 0, n 0 + 1, (71) a i (69) k {ˆx n,t }, t = 1,, k fundamental system (69) k t=1 d tx n,t (69) x n,t = rt n γ i rt i = 0. (7) r t P (r) = k γ ir i = 0 P (r) = 0 k r t {r n t } (69) fundamental system (68) x n = d t rt n + ψ n (7) t=1 P (r) = 0 r j, j = 1,, p µ j ( p j=1 µ j = k) x n = [d 1,1 + d 1, n + d 1, n(n 1) + + d 1,µ1 n(n 1) (n µ 1 1)]r n 1 + [d,1 + d, n + d, n(n 1) + + d,µ n(n 1) (n µ 1)]r n + + [d p,1 + d p, n + d p, n(n 1) + + d p,µ n(n 1) (n µ p 1)]r n p + ψ n (74) 8
trivial dx/dt = 0, x(0) = 0 α i x n i = 0 (75) ρ ζ ζ α i ζ n i = 0 (76) x n = t[d 1,1 + d 1, n + d 1, n(n 1) + + d 1,µ1 n(n 1) (n µ 1 1)]ζ1 n + t[d,1 + d, n + d, n(n 1) + + d,µ n(n 1) (n µ 1)]ζ n + + t[d p,1 + d p, n + d p, n(n 1) + + d p,µ n(n 1) (n µ p 1)]ζp n (77) t 0 x(t) = 0 lim t 0,n t=x tnq ζ n i = X lim n nq 1 ζ n i = 0 (78) ζ i < 1 Definition 1. ρ(ζ) 1 (zero stable) n e = 1 [ α i x(t n i ) t β i f(x(t n i ))] (79) α 0 = 1 α 0 [ α i x(t n i ) t β i x (t n i )] (80) m x m m p m p m (t) = (t t n ) m p 0,, p m m p j (j = 1,, m) e = 0 m (80) p j (t t n ) j = 0 α 0 e = α i = 0 (81) 9
j 0 α 0 e = = [α i (t n i t n ) j tβ i j(t n i t n ) j 1 ] (8) ( i t) j 1 [ i tα i j tβ i ] (8) = ( t) j [i j α i + ji j 1 β i ] = 0 (84) j=0 α i = 0 (85) i j α i = j i j 1 β i (j 1) (86) m α i, β i (i = 0,, k)k + 1 k + 1 m + 1 k m = k k > m = k + 1 [Dahlquist(1956)] m t n Taylor x(t) = x (t) = m j=0 p j (t)x (j) (t n ) j! m p j (t)x(j) (t n ) j=1 j! + p m+1(t)x (m+1) (ξ(t)) (m + 1)! + p m+1(t)x (m+1) (η(t)) (m + 1)! (87) (88) x (j) (t) x t j ξ(t),η(t) t n,t L e = Lx Lp 0 = Lp 1 = = Lp m = 0. [ ] p m+1 (t n i )x (m+1) (ξ i (t n i )) p Lx = α i tβ m+1(t n i )x (m+1) (η i (t n i )) i (89) (m + 1)! (m + 1)! i=1 i=1 t n i ξ i (t n i ), η i (t n i ) t n 1 p m+1 (t n i )x (m+1) (ξ i (t n i )) α i (m + 1)! α ( i t) m+1 x (m+1) (ξ i (t n i )) i (m + 1)! (90) i=1 tm+1 (m + 1)! α i i m+1 x (m+1) (ξ i (t n i)) (91) i=1 10
p β m+1(t n i )x (m+1) (η i (t n i )) i (m + 1)! β ( i t) m x (m+1) (η i (t n i )) i m! (9) i=1 i=1 tm m! β i i m x (m+1) (η i (t n i)) (9) i=1 Lx Lx = sup t [t n k,t n] x (m+1) (t) M t m+1 (94) M = 1 α i i m+1 + 1 β i i m (m + 1)! m! (95) i=1 i=1 Lx = O( t m+1 ) (96).1.1.1 1 k = 1 k = 1 k+1 = m+1 m = k = m m = 1 j = 1 α 0 + α 1 =0, (97) α 1 = (β 0 + β 1 ) (98) α 0 = 1 α 1 = 1, β 1 = 1 β 0 x n x n 1 = t (β 0 f n + (1 β 0 )f n 1 ). (99) () m = j = α 1 = β 1 (100) β 1 = 1/, β 0 = 1/ 1.1. Adams α 0 = 1 α i (i = 1,, k) 1 (α l ) ( α l = α 0 = 1) l = 1 Adams 1 Adams 1 β 0 11
Adams-Bashforth ( ( ) x n = x n 1 + t f(t n 1, x n 1 ) 1 ) f(t n, x n ), (101) ( ( ) x n = x n 1 + t 1 f(t n 1, x n 1 ) 4 f(t n, x n ) + 5 ) 1 f(t n, x n ). (10) k Adams-Bashforth α i =,..., k β β 0 = 0 k (k 1) 1 = k Adams-Moulton ( 5 ( ) x n = x n 1 + t 1 f(t n, x n ) + f(t n 1, x n 1 ) 1 ) 1 f(t n, x n ) (10) ( ( )4 x n = x n 1 + t 8 f(t n, x n ) + 19 4 f(t n 1, x n 1 ) (104) 5 4 f(t n, x n ) + 1 4 f(t n, x n ) Adams-Moulton β 0 k + 1 Nystöm: l = generalized Milen-Simpson: l = Cotes: k = l ) (105).1. BDF (Backward Differentiation Formulae) β i = 0 (i 0) Gear (BDF) (Stiffly-Stable Method) [Gear] β 0 = 1 α (85)(86) α i x n i = β 0 f(t n, x n ) (106) x n x n 1 + 1 x n = f(t n, x n ) (107) 11 6 x n x n 1 + x n 1 x n = f(t n, x n ). (108) k + 1 k k 6 1 BDF f(t, x) f [7]( BDF[explicit BDF] ) α i (i = 0,, k) β 0 (85)(86) β i (i = 1,, k ) k < i k i α i = 0 k = k k = 1 Adams-Bashforth 1
1: Coefficients of Gear s BDF Method Coefficient k = 1 k = k = k = 4 k = 5 k = 6 11 5 17 49 α 0 1 6 1 60 0 α 1 1 4 5 6 1 15 α 0 5 α 0 0 1 4 10 0 α 4 0 0 0 1 4 α 5 0 0 0 0 1 5 6 5 α 6 0 0 0 0 0 1 6 5 4 15 4 Coefficient k = 1 k = 1 : Coefficients of explicit BDF Method k = k = k = k = k = k = 4 k = k = k = 4 k = 4 11 5 α 0 1 6 1 α 1 1 4 1 1 1 α 0 α 0 0 0 0 1 4 1 α 4 0 0 0 0 0 4 β 0 0 0 0 0 0 0 8 1 β 1 1 4 4 β 0 1 7 49 β 0 0 1 6 9 1 1 4 β 4 0 0 0 7 1 0 1 1
.1.4 (Predictor-Corrector) 1 1) (P) )f(x n ) (E) ) (C) 4) x n f(x n ) (E) ᾱ i x n i = t α i x n i = t β i f(t n i, x n i ) (109) i=1 β i f(t n i, x n i ) (110) 4 PEC α 0 x (0) n + ᾱ i x i = t β i f(t i, x (0) i ) (111) i=1 α i x i = t i=1 β i f(t i, x (0) i ) (11) 1 PEC(E) P(EC) (E).1.5 Lagrange Polynomial Interpolation x f N +1 (x j, f(x j )), j = 0,, N f(x) N P N P N (x) = l 0 0(x) =1, N f(x j )l N j (x) (11) j=0 l N j (x) = (x x 0)(x x 1 ) (x x j 1 )(x x j+1 ) (x x N ) (x j x 0 )(x j x 1 ) (x j x j 1 )(x j x j+1 ) (x j x N ) (114) (N 1). (115) (t n i, f(t n i, x n i )) i = 1,..., N P e N (t n, f(t n, x n )) e explicit (t n 1, f(t n 1, x n 1 )) P1 e f x n x n 1 + tn t n 1 P e 1 (t)dt = x n 1 + tn t n 1 f(t n 1, x n 1 )dt = x n 1 + (t n t n 1 )f(t n 1, x n 1 ) (116) Euler P e, P e Adams-Bashforth P e t t n t t n 1 (t) =f(t n 1, x n 1 ) + f(t n, x n ) (117) t n 1 t n t n t n 1 P e (t t n )(t t n ) (t) =f(t n 1, x n 1 ) (t n 1 t n )(t n 1 t n ) + f(t (t t n 1 )(t t n ) n, x n ) (t n t n 1 )(t n t n ) (t t n 1 )(t t n ) + f(t n, x n ) (t n t n 1 )(t n t n ). (118) 14
tn x n x n 1 P e (t)dt t n 1 tn ( ) t t n t t n 1 = f(t n 1, x n 1 ) + f(t n, x n ) dt t n 1 t n 1 t n t n t n 1 [ ( ) f(tn 1, x n 1 ) t = t n 1 t n t n t f(t ( )] n, x n ) t tn t n 1 t n t n 1t t n 1 = f(t ( ) n 1, x n 1 ) t0 (t n + t n 1 ) t n t 0 f(t ( ) n, x n 1 ) t0 (t n + t n 1 ) t n 1 t 0 t 1 t 1 =f(t n 1, x n 1 ) t 0 t 1 ((t n t n 1 ) + (t n 1 t n )) f(t n, x n ) t 0 t 1 ((t n t n 1 )) =f(t n 1, x n 1 ) t 0 ( t 0 + t 1 ) t 1 f(t n, x n ) t 0 t 1. (119) t i = t n i t n i 1 i t = t i ( x n x n 1 t f(t n 1, x n 1 ) 1 ) f(t n, x n ) (10) 15
tn x n x n 1 P e (t)dt t n 1 tn ( (t t n )(t t n ) = f(t n 1, x n 1 ) t n 1 (t n 1 t n )(t n 1 t n ) + f(t (t t n 1 )(t t n ) n, x n ) (t n t n 1 )(t n t n ) ) (t t n 1 )(t t n ) +f(t n, x n ) dt (t n t n 1 )(t n t n ) [ ( )] f(t n 1, x n 1 ) t tn = (t n 1 t n )(t n 1 t n ) (t n + t n ) t + t n t n t t n 1 [ ( )] f(t n, x n ) t tn + (t n t n 1 )(t n t n ) (t n 1 + t n ) t + t n 1t n t t n 1 [ ( )] f(t n, x n ) t tn + (t n t n 1 )(t n t n ) (t n 1 + t n ) t + t n 1t n t t n 1 ( = f(t ( ) n 1, x n 1 ) t0 t n + t n t n 1 + tn 1) (t n + t n ) t 0 (t n + t n 1 ) + t n t n t 0 t 1 ( t 1 + t ) ( + f(t ( ) n, x n ) t0 t n + t n t n 1 + tn 1) (t n 1 + t n ) t 0 (t n + t n 1 ) + t n 1 t n t 0 t 1 t ( + f(t ( ) n, x n ) t0 t n + t n t n 1 + tn 1) (t n 1 + t n ) t 0 (t n + t n 1 ) + t n 1 t n t 0 ( t 1 + t ) t = f(t n 1, x n 1 ) t 0 ( ( t 6 t 1 ( t 1 + t ) n + t n t n 1 + t ) n 1) (tn + t n ) (t n + t n 1 ) + 6t n t n f(t n, x n ) t 0 6 t 1 t + f(t n, x n ) t 0 6 t ( t 1 + t ) ( ( t n + t n t n 1 + t ) ) n 1 (tn 1 + t n ) (t n + t n 1 ) + 6t n 1 t n ( ( t n + t n t n 1 + t ) ) n 1 (tn 1 + t n ) (t n + t n 1 ) + 6t n 1 t n. (11) t i = t n t n i = i 1 j=0 t j (t n + t n t n 1 + t n 1) (t n i + t n j ) (t n + t n 1 ) + 6t n i t n j ( = t n + t n (t n t 0 ) + (t n t 0 ) ) ( t n t i t j) (tn t 0 ) + 6 (t n t i) ( t n t ) j = t 0 t 0 ( t i + t j) + 6 t i t j (1) 16
t n x n x n 1 t ( 0 t 0 t 0 ( t 0 + t 1 + t ) + 6 ( t 0 + t 1 ) ( t 0 + t 1 + t ) ) f(t n 1, x n 1 ) 6 t 1 ( t 1 + t ) t ( 0 t 0 t 0 ( t 0 + t 1 + t ) + 6 t 0 ( t 0 + t 1 + t ) ) f(t n, x n ) 6 t 1 ( t 1 + t ) + t ( 0 t 0 t 0 ( t 0 + t 1 ) + 6 t 0 ( t 0 + t 1 ) ) f(t n, x n ) 6 t ( t 1 + t ) = t 0 + t 0 ( t 1 + t ) + 6 t 0 t 1 ( t 1 + t ) f(t n 1, x n 1 ) 6 t 1 ( t 1 + t ) t 0 + t 0 ( t 1 + t ) 6 t 1 t f(t n, x n ) + t 0 + t 0 t 1 6 t ( t 1 + t ) f(t n, x n ) (1) x n x n 1 1 f(t n 1.x n 1 ) 4 f(t n, x n ) + 5 1 f(t n, x n ) (14) Adams-Bashforth (t n, f(t n, x n )) N + 1 P i N P i 0 =f(t n, x n ), (15) P1 i =f(t n, x n ) t t n 1 t t n + f(t n 1 ), (16) t n t n 1 t n 1 t n P i =f(t n, x n ) (t t n 1)(t t n ) (t n t n 1 )(t n t n ) + f(t (t t n )(t t n ) n 1, x n 1 ) (t n 1 t n )(t n 1 t n ) (t t n )(t t n 1 ) + f(t n, x n ) (17) (t n t n )(t n t n 1 ) x n x n 1 t 0 f(t n, x n ), (18) x n x n 1 t 0 (f(t n, x n ) + f(t n 1, x n 1 )), (19) x n x n 1 0 + t 0 t 1 6( t 0 + t 1 ) f(t n, x n ) + t 0 + t 0 t 1 f(t n 1, x n 1 ) 6 t 1 Euler Adams-Moulton t 0 6 ( t 0 + t 1 ) t 1 f(t n, x n ). f(t n, x n ) x BDF (t n i, x n i ) (i = 0,..., N) x Q N (t) (10) t t n 1 t t n Q 1 (t) = x n + x n 1 (11) t n t n 1 t n 1 t n dx dt dq 1 = dt x n t n t n 1 + x n 1 t n 1 t n = x n x n 1 t 0 = f(t n, x n ). (1) 17
(t t n 1 )(t t n ) Q (t) = x n (t n t n 1 )(t n t n ) + x (t t n )(t t n ) n 1 (t n 1 t n )(t n 1 t n ) + x (t t n )(t t n 1 ) n (t n t n )(t n t n 1 ) dx dt dq dt t (t n 1 + t n ) =x n (t n t n 1 )(t n t n ) + x t (t n + t n ) n 1 (t n 1 t n )(t n 1 t n ) + x t (t n + t n 1 ) n (t n t n )(t n t n 1 ). (14) t = t n dx t n t n 1 t n dt x n t=tn (t n t n 1 )(t n t n ) + x t n t n t n n 1 (t n 1 t n )(t n 1 t n ) + x t n t n t n 1 n (t n t n )(t n t n 1 ) t 0 = t 1 (1) t 0 + t 1 =x n t 0 ( t 0 + t 1 ) x t 0 + t 1 t 0 n 1 + x n = f(t n, x n ). (15) t 0 t 1 ( t 0 + t 1 ) t 1 4 Q x n x n 1 + 1 x n = f(t n, x n ). (16) (t t n 1 )(t t n )(t t n ) Q (t) =x n (t n t n 1 )(t n t n )(t n t n ) + x (t t n )(t t n )(t t n ) n 1 (t n 1 t n )(t n 1 t n )(t n 1 t n ) (t t n )(t t n 1 )(t t n ) + x n (t n t n )(t n t n 1 )(t n t n ) + x (t t n )(t t n 1 )(t t n ) n (t n t n )(t n t n 1 )(t n t n ). (17) x Q t = t n dx x n ( dt t t=tn t 0 ( t 0 + t 1 )( t 0 + t 1 + t ) n t n (t n 1 + t n + t n ) + (t n 1 t n + t n t n + t n t n 1 ) ) x n 1 ( t t 0 t 1 ( t 1 + t ) n t n (t n + t n + t n ) + (t n t n + t n t n + t n t n ) ) x n ( + t ( t 0 + t 1 ) t 1 t n t n (t n + t n 1 + t n ) + (t n t n 1 + t n 1 t n + t n t n ) ) x n ( t ( t 0 + t 1 + t )( t 1 + t ) t n t n (t n + t n 1 + t n ) + (t n t n 1 + t n 1 t n + t n t n ) ) x n = t 0 ( t 0 + t 1 )( t 0 + t 1 + t ) ( t 0( t 0 + t 1 ) + ( t 0 + t 1 )( t 0 + t 1 + t ) + t 0 ( t 0 + t 1 + t )) x n 1 t 0 t 1 ( t 1 + t ) ( t 0 + t 1 )( t 0 + t 1 + t ) x n + t 0 ( t 0 + t 1 + t ) ( t 0 + t 1 ) t 1 t x n t 0 ( t 0 + t 1 ) ( t 0 + t 1 + t )( t 1 + t ) t =x n t 0 + (4 t 1 + t ) t 0 + t 1 ( t 1 + t ) t 0 ( t 0 + t 1 )( t 0 + t 1 + t ) x n 1 ( t 0 + t 1 )( t 0 + t 1 + t ) t 0 t 1 ( t 1 + t ) + x n t 0 ( t 0 + t 1 + t ) ( t 0 + t 1 ) t 1 t x n t 0 ( t 0 + t 1 ) ( t 0 + t 1 + t )( t 1 + t ) t =f(t n, x n ). (18) 18
t 0 = t 1 = t 11 6 x n x n 1 + x n 1 x n = f(t n, x n ) (19) dq dt t 0 + t 1 x n t 0 ( t 0 + t 1 ) x t 0 + t 1 t 0 n 1 + x n t 0 t 1 ( t 0 + t 1 ) t 1 P e (t n ) (140) tn = f(t n 1, x n 1 ) t n t n t n 1 t n + f(t n, x n ) t n t n 1 t n t n 1 dq dt x n t 0 + (4 t 1 + t ) t 0 + t 1 ( t 1 + t ) t 0 ( t 0 + t 1 ) ( t 0 + t 1 + t ) = f(t n 1, x n 1 ) t 0 + t 1 t 1 f(t n, x n ) t 0 t 1, (141) P e (t n ) (14) t=tn x n 1 ( t 0 + t 1 ) ( t 0 + t 1 + t ) t 0 t 1 ( t 1 + t ) t 0 ( t 0 + t 1 + t ) t 0 ( t 0 + t 1 ) + x n x n ( t 0 + t 1 ) t 1 t ( t 0 + t 1 + t ) ( t 1 + t ) t (t n t n )(t n t n ) = f(t n 1, x n 1 ) (t n 1 t n )(t n 1 t n ) + f(t (t n t n 1 )(t n t n ) n, x n ) (t n t n 1 )(t n t n ) (t n t n 1 )(t n t n ) + f(t n, x n ) (t n t n 1 )(t n t n ) = f(t n 1, x n 1 ) ( t 0 + t 1 )( t 0 + t 1 + t ) t 1 ( t 1 + t ) t 0 = t 1 = t f(t n, x n ) t 0( t 0 + t 1 + t ) t 1 t + f(t n, x n ) t 0( t 0 + t 1 ) ( t 1 + t ) t (14) 11 6 x n x n 1 + x n 1 x n = f(t n 1, x n 1 ) f(t n, x n ) + f(t n, x n ). (144) dq dt t 0 + t 1 x n t 0 ( t 0 + t 1 ) x t 0 + t 1 t 0 n 1 + x n t 0 t 1 ( t 0 + t 1 ) t 1 = f(t n 1, x n 1 ) ( t 0 + t 1 )( t 0 + t 1 + t ) t 1 ( t 1 + t ) P e (t n ) (145) t=tn f(t n, x n ) t 0( t 0 + t 1 + t ) t 1 t + f(t n, x n ) t 0( t 0 + t 1 ) ( t 1 + t ) t (146) 19
(85), (86) dq dt c 1 P e (t n ) + c P e (t n ) (147) t=tn (c 1 + c = 1) t 0 + t 1 x n t 0 ( t 0 + t 1 ) x t 0 + t 1 t 0 n 1 + x n t 0 t 1 ( t 0 + t 1 ) t 1 = f(t n 1, x n 1 ) ( t 0 + t 1 )(c t 0 + t 1 + t ) t 1 ( t 1 + t ) c = / t 0 = t 1 = t f(t n, x n ) t 0(c ( t 0 + t 1 ) + t ) t 1 t + f(t n, x n ) c t 0 ( t 0 + t 1 ) ( t 1 + t ) t. (148) x n x n 1 + 1 x n = 8 f(t n 1, x n 1 ) 7 f(t n, x n ) f(t n, x n ) (149) c = /. Taylor Hamilton ( ) Symplectic(leap-flog) 4 bound bound ( (standard test problem) ) dx dt = λx, λ C. (150) x Rλ 0 bound (Rλ < 0 ) Euler x n = x 0 (1 + λ t) n λ t 1 [bounded] ( λ t < 1 [convergent]) λ t bound λ t (stability region) ( strict stability region) Definition (A-stability (Dahlquist)). C {λ : Rλ < 0} A 0
A Runge-Kutta A A A A R R Runge-Kutta A A C(α) = {λ C : arg(λ) α} A(α) [Widlund] (Stiff Stability)[Gear] Definition (Widlund). α (0, π/) W α = {λ t : π arg(λ t) < α} A(α) α A(0) A(0) k k + 1 α [0, π/] k = m =, k = m = 4 k m A(α) Definition 4 (Gear). R 1, R [stiffly stable] R 1 = {λ t : Rλ t < a} (151) R = {λ t : a Rλ t b, Iλ t c} (15) a, b, c λ t R 1 4.1 Runge-Kutta x 0 = 1 Runge-Kutta x n = ξ(λ t) n K = (k (1), k (),, k (R) ) T, B = (b rs ), C = (c r ), E = (1, 1, ) K = λ(x n 1 E + tbk) (15) x n = x n 1 + tc T K (154) K x n = x n 1 + tc T (I λ tb) 1 λx n 1 E. (155) ξ(λ t) = 1 + λ tc T (I λ tb) 1 E. (156) ξ(λ t) 1 4.1 Runge-Kutta (R 4) R > 4 6 5 Runge-Kutta Lawson[9] 1
Im(λ t) R1 c R -a b Re(λ t) -c 4. (150) (α i λ tβ i )x n i = 0. (157) ρ(ζ) = σ(ζ) = α i ζ i (158) β i ζ i (159) π(ζ, λ t) = ρ(ζ) λ tσ(ζ) = 0 (160) ζ i (157) x n = t[d 1,1 + d 1, n + d 1, n(n 1) + + d 1,µ1 n(n 1) (n µ 1 1)]ζ1 n + t[d,1 + d, n + d, n(n 1) + + d,µ n(n 1) (n µ 1)]ζ n + + t[d p,1 + d p, n + d p, n(n 1) + + d p,µ n(n 1) (n µ p 1)]ζp n (161) ζ i < 4. 4.
4 Stability Diagram of Runge-Kutta Method Stage[ R]& Order[ m ] R=1, m=1 R=, m= R=, m= R=4, m=4 R=6, m=5[ Lawson] Im( λ t) 0 - -4-7 -6-5 -4 - - -1 0 Re( λ t) 1: Stability diagram of Runge-Kutta method.
1.5 1 Stability Diagram of Adams-Bashforth Method Number of Steps k=1 k= k= k=4 k=5 k=6 0.5 Im( λ t) 0-0.5-1 -1.5-1.5-1 -0.5 0 0.5 1 1.5 Re( λ t) : Stability diagram of Adams-Bashfortht method. 4
4 Stability Diagram of Adams-Moulton Method Number of Steps k=1 k= k= k=4 k=5 k=6 1 Im( λ t) 0-1 - - -6-5 -4 - - -1 0 1 Re( λ t) : Stability diagram of Adams-Moulton method. 5
0 15 10 Stability Diagram of Gear s BDF Method Number of Steps k=1 k= k= k=4 k=5 k=6 5 Im( λ t) 0-5 -10-15 -5 0 5 10 15 0 5 0 Re( λ t) 4: Stability diagram of BDF method. 6
1.5 1 Stability Diagram of explicit BDF Method Number of Steps k=1 k= k= k=4 k=[ α k=[ =0] =0, α 4 =0] 0.5 Im( λ t) 0-0.5-1 - -1.5-1 -0.5 0 0.5 Re( λ t) 5: Stability diagram of explicit BDF method. 7
4. stiffness Definition 5. λ t, t = 1,,, m dx dt = Ax + Φ(t) (16) Rλ t < 0, t = 1,,, m max t=1,,,m Rλ t min t=1,,,m Rλ t (stiff) stiffness ratio [1] J.C. Butcher, The Numerical Analysis of Ordinary Differential Equations (John Wiley & Sons, Chichester/New York, 1987). [] J.D. Lambert, Computational Methods in Ordinary Differential Equations (John Wiley & Sons, London, 197). [] C.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations (Prentice-Hall, Englewood Cliffs, N.J., 197). [4] R.D. Richtmyer, K.W. Morton, Difference Methods for Initial Value Problems (Interscience, New York, 1957). [5] P. Hartman, Ordinary Differential Equations (John Wiley & Sons, New York, 1964). [6] C.A.J. Fletcher, Computational Techniques for Fluid Dynamics (Springer, Berlin, 1991). [7] G.E. Karniadakis et al, J. Comput. Phys 97, 414 (1991); M.A. Hulsen, MEAH-18 (1996). [8] S. Gill, Proc. Cambridge Philos. Soc. 47, 95 (1951). [9] J.D. Lawson, SIAM J. Numer. Anal. 9??(197). 8