17 1 25

http://grape.astron.s.u-tokyo.ac.jp/pub/people/makino/kougi/stellar_dynamics/index.html http://grape.astron.s.u-tokyo.ac.jp/pub/people/makino/talks/index-j.html

d 2 x i dt 2 = j i Gm j x j x i x j x i 3, (1) 2 10

2 O(N 2 ) N

( )

k 1 = x i + h 2 f(x i, t i ) x i+1 = x i + hf(k 1, t i + h/2) (2) y step 2 y i step 1 h x i x x i+1 2 ( 2 )

x n+1 k i = x n + h s i=1 = f(x n + h s b i k i j=1 a ij k j, t n + c i h) (3) s : number of stages a ij, b i, c i a c c i = s j=1 a ij (4)

2 s = 2 x n+1 = x n + h(k 1 b 1 + k 2 b 2 ) k 1 = f(x n + h(a 11 k 1 + a 12 k 2 ), t n + c 1 h) k 2 = f(x n + h(a 21 k 1 + a 22 k 2 ), t n + c 2 h) (5) a ij (j i) 0 k 1 a ij (j > i) 0 k i k i (semi-implicit) k 1 k 2

k i

Runge-Kutta x n+1 = x n + h(k 1 /6 + k 2 /3 + k 3 /3 + k 4 /6) k 1 = f(x n, t n ) k 2 k 3 = f(x n + hk 1 /2, t n + h/2) = f(x n + hk 2 /2, t n + h/2) k 4 = f(x n + hk 3, t n + h) (6)

Runge-Kutta 1. a ij i j = 1 0 2. 4 4 3. 4

s 1 2 3 4 5 6 7 8 9 10 s 10 p 1 2 3 4 4 5 6 6 7 7 p s 2 2 2

1. 2. 3. 0 Dormand Prince http://www.unige.ch/math/folks/hairer/software.html (Fotran C )

: : (variable step) / (adaptive step) h =

RK x n+1 = x n + s i=1 b i k i (7) k i b i b i x = x n + s i=1 b i k i (8) x embedded Runge-Kutta-Fehlberg Dormand-Prince

( ) s 2s 8 2s

: t i x i t i+1 x i+1 : t i

: f f i-3 i-2 i-1 i i+1 t

i i + 1 i p i p P (t j ) = f j = f(x j, t j ) (i p j i) (9) i + 1 x i+1 x i+1 = x i + t i+1 t i P (t)dt (10) h p x i+1 = x i + h p l=0 a pj f i l (11)

2 p = 1 P (t) = f i f i 1 f i (t t i ) (12) h x i+1 = x i + h 2 (3f i f i 1 ) (13)

t i+1 h p x i+1 = x i + h p 1 l= 1 b pj f i l (14) p = 1 x i+1 = x i + h 2 (f i + f i+1 ) (15)

- - (PEC) - (PE) P(EC) 2 PE(CE) 2

t 0 x 0 : : :

: :

2

( )

2 2 : 2

: : ( )

: ρ r α σ r (2 α)/2 (16) : 1/ Gρ r α/2

2 0 2 : : 0

: : O(N) 2 : O(N 1/3 )

( ) 2

t i t i 1. t i + t i 1 2. 2 3. 4. 1 i t i ti : t i + t i n Time

PEC

t = η a 1 a (2) 1 + ȧ 1 2 ȧ 1 a (3) 1 + a (2) (17) 1 2

:

3 :? /

: 2 leap frog leapfrog Verlet 2nd order ABM 2nd order Stömer-Cowell 4

leap frog v i+1/2 = v i 1/2 + ta(x i ) (18) x i+1 = x i + tv i+1/2 (19) ±1/2 v 1/2 = v + ta(x 0 )/2 (20) v i = v i 1/2 + ta(x i )/2 (21)

x i+1 = x i + tv i + t 2 a(x i )/2 (22) v i+1 = v i + t[a(x i ) + a(x i+1 )]/2 (23) PC v i+1/2 = v i + ta(x i )/2 (24) x i+1 = x i + tv i+1/2 (25) v i+1 = v i+1/2 + ta(x i+1 )/2 (26) x i 1, x i, x i+1

leap frog 2 :

d 2 x = x (27) dt2 2 (x, v) = (1, 0) 1/4

H = 1 2 (x2 + v 2 ) h2 4 x2 (28) h 2 ( 1 )

d 2 x dt = 2 x3 (29) 2 (x, v) = (1, 0) 1/8

1. 2.

1. 2. 3.

1. 2.

? : :

: 7 6 15 8 T (p) + V (q)

Ruth 4 S 4 (h) = L(d 1 h)l(d 2 h)l(d 1 h), d 1 = 1/(2 2 1/3 ), d 2 = 1 2d 1 = 2 1/3 /(2 2 1/3 ) (30) : L(h) h 1

4 h t t i i+1 2 1/3 2 ( ) ( )

6 6 9 8 27 6 : d 1 = d 7 = 0.784513610477560 d 2 = d 6 = 0.235573213359357 d 3 = d 5 = 1.17767998417887 d 4 = 1.31518632068391 (31)

? H p H H H H = H + h p+1 H p + O(h p+2 ) (32) H ( ) ( )

( ) H H = H 1 + H 2 + (33)

1 : ( ) :8 15 2 Dormand 9 8 8 10 4 100

MVS + ( ) : H = T (p) + V (q) (34) p p p t H q (35) q

: H = T (p) + V 1 (q) + V 2 (q) (36) V 1 V 2 H 1 = T (p) + V 1 (q) (37) V 2 p p t V 2 (38) q 2

: x x + tv MVS: V 1 KS

( ) dy dx = f(x, y) (39) y i+1 = F (x i, y i, f, x) (40) y i = F (x i+1, y i+1, f, x) (41) 1

y i+1 = y i + h 2 (f i + f i+1 ) (42)

:

RK

= 2 : α 0 x i+p + α p x i = h 2 (β 0 f i+p + β p f i ) (43) β 0 = 0 α i = α p i, β i = β p i (44) x i+1 2x i + x i 1 = h 2 f i (45)

(x i,...x i+p 1 ) x i+p (x i+p,...x i+1 ) x i

1970 1990 Quinlan Tremaine 14 Quinlan & Tremaine

: 6 6

: 4 x 1 = x 0 + t 2 (v 1 + v 0 ) t2 12 (a 1 a 0 ), (46) v 1 = v 0 + t 2 (a 1 + a 0 ) t2 12 (ȧ 1 ȧ 0 ). (47) ( )

( )

e = 0.9

t 0 t 0 t 1 t 1 t = 2 1[h(ξ 0) + h(ξ 1 )] (48) h(ξ) ξ 1 t

Cano Sanz-Serna (a), (b) (c) 1