Accurate Numerical Integration Algorithm for the Kepler Motion

Accurate Numerical Integration Algorithm
for the Kepler Motion
Yosuke NAKANISHI (Pixela Corporation),
Yoshimasa NAKAMURA (Kyoto University),
Yukitaka MINESAKI (Tokushima Bunri University)
Annual Workshop on Modeling and Simulation in Applied Mathematics in Josai University
Tokyo, Japan, December. 20. 2015
Outline
1
Kepler Motion
1
Existing Numerical Integrators and Kepler Motion
2
Discrete Gradient Method
3
Kepler Motion and Harmonic Oscillator
Totally Conservative Integrator for Kepler Motion (Discrete-time Kepler Motion)
4
1
5
6
Reproducibility of Orbits in Discrete Kepler Motion
Time Adjustment of Discrete Kepler Motion
Conclusion
1
Kepler 運動
互い万有引力を及ぼし合う 2 質点間の相互運動
運動方程式
dx
dpx
x
= px ,
= −K 2 3 ,
dt
dt
|x|
!
K = G(m1 + m2 ),
m1 , m2 : 質点 1,2 の質量,
G : 万有引力定数
保存量
1
K2
2
エネルギー = Hamiltonian :
E = H = |px | −
2
|x|
角運動量 :
h = x × px
x
Runge-Lenz (Laplace) ベクトル : e = px × h − K 2
|x|
2
数値積分法による Kepler 運動の再現
数値積分法がもつ保存性
symplectic
エネルギー
全保存型差分法
数値積分法
保存差分法
(YM & YN, 2002)
×
⃝
×
⃝
×
×
⃝
⃝
⃝
エネルギー
角運動量
Runge-Lenz ベクトル
z
E=-0.5
E=-0.2
E= 0.0
E= 0.2
E= 0.5
z
1
0.5
0
-0.5
0
-1
-1
-1.5
-2
-2
-2.5
-3
0.5
-4-1.5
-0.5
-1
-0.5
x
-1
0
0.5
-3
0.5
-3.5
0
-4-3.5
y
0
-3
-0.5
-2.5
-2
-1.5
-1.5
x
1 -2
symplectic 数値積分法による軌道の再現
-1
-1
-0.5
0
-1.5
0.5
1 -2
全保存型差分法よる軌道の再現
3
y
離散変分法 (1)
定義 (離散変分)
(j,j+1)
離散変分 ∇p
H(p, q), ∇q
は以下の関係式を満たす．
(i) (p(j+1)−p(j) ) · ∇p
(j,j+1)
(j,j+1)
+(q (j+1)−q (j) ) · ∇q
H(p, q)
H(p, q)
(j,j+1)
H(p, q)
= H(p(j+1) , q (j+1) ) − H(p(j) , q (j) )
(ii) ∇p
∇q
(j,j)
(j,j)
H(p, q) = ∇p H(p(j) , q (j) ),
H(p, q) = ∇q H(p(j) , q (j) )
p = (p1 , · · · , pn ), q = (q1 , · · · , qn );
Hamiltonian (= エネルギー) : H(p, q);
•(j) は時刻 t(j) における • の値を示す．
4
離散変分法 (2)
本発表では，以下の Gozalez による離散変分 (1996) を用いる．
∇p
(j,j+1)
H(p, q)
!
"
(j+1)
(j)
(j+1)
(j)
p
+p
q
+q
H(p(j+1),q (j+1) )−H(p(j),q (j) )
:= ∇p H
,
+ (j+1)
(p(j+1) −p(j) )
(j)
(j+1)
(j)
2
2
2
2
|p
+p | + |q
+q |
#
$
(j+1)
(j)
(j+1)
(j)
p
+
p
q
+
q
(p(j+1) −p(j) ) · ∇p H
,
2
2
−
(p(j+1) −p(j) )
(j+1)
(j)
(j+1)
(j)
2
2
|p
+p | + |q
+q |
∇q
(j,j+1)
H(p, q)
!
"
(j+1)
(j)
(j+1)
(j)
p
+p
q
+q
H(p(j+1),q (j+1) )−H(p(j),q (j) )
:= ∇q H
,
+ (j+1)
(q (j+1) −q (j) )
(j)
(j+1)
(j)
2
2
2
2
|p
+p | + |q
+q |
#
$
p(j+1) +p(j) q (j+1) +q (j)
(q (j+1) −q (j) ) · ∇q H
,
2
2
−
(q (j+1) −q (j) )
(j+1)
(j)
(j+1)
(j)
2
2
|p
+p | + |q
+q |
5
離散変分法 (3)
1
エネルギー保存型数値積分法
1
t(n+1) −t(n)
→0
⎧
⎪
q (j+1) −q (j)
=⇒
(j,j+1)
⎪
⎨
=
∇
H(p,
q)
p
t(n+1) −t(n)
p(j+1) −p(j)
⎪
(j,j+1)
⎪
⎩ (n+1) (n) = −∇q
H(p, q)
t
−t
=⇒ エネルギー (= Hamiltonian) H(p, q) が保存される．
(代入すると，離散変分の定義 (i) の左辺が 0 になる．)
正準方程式
⎧
⎪
⎨ dq = ∇p H(p, q)
dt
dp
⎪
⎩
= −∇q H(p, q)
dt
しかしながら，一般にエネルギー以外は保存されない．
=⇒ 厳密には，2 以上の自由度をもつ Hamilton 系の軌道を再現できない．
エネルギー
角運動量
Runge-Lenz ベクトル
symplectic
エネルギー
全保存型差分法
数値積分法
保存差分法
(YM & YN, 2002,2004)
×
⃝
×
⃝
×
×
⃝
⃝
⃝
6
2 次元 Kepler 運動と 2 重調和振動子 (1)
2 重調和振動子 =⇒ 2 次元 Kepler 問題 (2D-Kepler)
dQk
dPk
2 つの (1 次元) 調和振動子 :
= Pk ,
= 2h2d Qk , k = 1, 2
ds
ds
1
エネルギーの値 : 2K 2 = (P12 + P22 ) − h2d (Q21 + Q22 ), 時刻 : s
2
⇓ 時間変換 : s → t :
ds
1
= 2
, k = 1, 2
dt
Q1 + Q22
LC 座標系上の 2D-Kepler :
dQk
Pk
dPk
2h2d Qk
= 2
,
=
, k = 1, 2,
dt
Q1 + Q22 dt
Q21 + Q22
1 P12 + P22
2K 2
エネルギーの値 : h2d =
− 2
, 時刻 : t
2
2
2
2 Q1 + Q2
Q1 + Q2
'
1& 2
2
⇓ Levi-Civita 変換 : x = Q1 Q2 , y =
Q1 − Q2 ,
2
P1 Q2 + P2 Q1
P1 Q1 − P2 Q2
px =
,
p
=
y
Q21 + Q22
Q21 + Q22
7
2 次元 Kepler 運動と 2 重調和振動子 (2)
2 重調和振動子 =⇒ 2 次元 Kepler 問題 (2D-Kepler)
相対座標系上の 2D-Kepler :
dx
dpx
x
= px ,
= −K 2
3 ,
2
2
dt
dt
(x + y ) 2
dy
dpy
y
= py ,
= −K 2
3 ,
2
2
dt
dt
(x + y ) 2
'
1& 2
K2
2
エネルギーの値 : h2d =
p + py − !
, 時刻 : t
2
2
2 x
x +y
まとめ時間変換，正準変換によって，
“2 つの 1 次元の振動子からなる系” は 2D-Kepler に変換される．
8
2 次元 Kepler 運動と 2 重調和振動子 (3)
1
保存量の関係
2 重調和振動子
相空間の次元 = 4 (自由度 2)
1
独立な保存量 (3 個)
1
2
3
調和振動子 1 のエネルギー
P12
H1 :=
− h2d Q21
2
調和振動子 2 のエネルギー
P22
H2 :=
− h2d Q22
2
ただし，2K 2 = H1 + H2
角運動量
LC 座標系上の 2D-Kepler
相空間の次元 = 4 (自由度 2)
1
独立な保存量 (3 個)
1
エネルギー
1 P12 + P22
2K 2
E = H :=
− 2
2 Q21 + Q22
Q1 + Q22
2
RL ベクトルの y 成分
(P1 Q2 −P2 Q1 )(P1 Q2 +P2 Q1 )+2K 2 (Q21 −Q22 )
ey := −
2(Q21 + Q22 )
3
角運動量の z 成分
hz := P1 Q2 − P2 Q1
J := P1 Q2 − P2 Q1
RL ベクトルの x 成分は E, ey , hz の関数
9
2 次元 Kepler 運動と 2 重調和振動子 (4)
1
保存量の関係
2 重調和振動子
相空間の次元 = 4 (自由度 2)
1
独立な保存量 (3 個)
1
2
調和振動子 1 のエネルギー
P12
H1 :=
− h2d Q21
2
調和振動子 2 のエネルギー
P22
H2 :=
− h2d Q22
2
LC 座標系上の 2D-Kepler
相空間の次元 = 4 (自由度 2)
1
独立な保存量 (3 個)
1
エネルギー
E = H :=
2
H1 + H2 − 2K 2
Q21 + Q22
+ h2d = h2d
RL ベクトルの y 成分
1
1
ey := − H1 + H2
2
2
ただし， 2K 2 = H1 +H2
3
角運動量
3
角運動量の z 成分
J := P1 Q2 − P2 Q1
hz := J
まとめ 2 重調和振動子がもつ独立の保存量 H1 , H2 , J が保たれるならば，
2D-Kepler がもつ全ての保存量も保たれる．
10
2 次元 Kepler 運動の全保存型差分 (1)
調和振動子 1,2 のエネルギー H1 , H2 ; 時刻 s
⇓ 離散変分法 (3)
離散 2 重調和振動子
(j+1)
(j)
(j+1)
(j)
(j+1)
(j)
(
)
Qk
−Qk
Pk
+Pk
Pk
−Pk
(j+1)
(j)
=
, (j+1)
= h2d Qk
+Qk ,
2
s(j+1) − s(j)
s
− s(j)
k = 1, 2; j = 0, 1, · · · ; 離散時刻 : s(n) , s(n+1)
1
独立な保存量 : H1 , H2 , J の 3 個
1
⇓ 離散時間変換 s
1
(j)
→t
(j)
s(j+1) − s(j)
1
: (j+1)
=
(j)
(j)
t
− t(j)
(Q1 )2 + (Q2 )2
離散 2D-Kepler
(j+1)
(j)
(j+1)
(j)
(j+1)
(j)
(j+1)
(j)
Qk
−Qk
+Pk
Pk
−Pk
Qk
+Qk
1 Pk
=
,
=
h
,
2d
(j) 2
(j+1) − t(j)
(j) 2
(j) 2
2
t(j+1) − t(j) 2 (Q(j)
t
(Q1 ) +(Q2 )
1 ) +(Q2 )
k = 1, 2; j = 0, 1, · · · ; 離散時刻 : t(n) , t(n+1)
1
独立な保存量 : H, ey (RL ベクトルの y 成分), hz (角運動量の z 成分) の 3 個
11
2 次元 Kepler 運動の全保存型差分 (2)
離散 2D-Kepler において，
(運動の自由度)
1
=
=
(相空間の次元)
4
−
−
(独立な保存量の数)
3
が成り立つが，独立な保存量は全て，元の 2D-Kepler と同じもの．
=⇒ 離散 2D-Kepler の軌道は 2D-Kepler の軌道に一致する．
h2d < 0 (ϵ < 1) : 楕円軌道
h2d = 0 (ϵ = 1) : 放物線軌道
h2d > 0 (ϵ > 1) : 双曲線線軌道
2 次元 Kepler 運動の全保存型差分
12
3 次元 Kepler 運動と 4 重調和振動子 (1)
4 重調和振動子 + 1 拘束条件 ⇒ 3 次元 Kepler 問題 (3D-Kepler)
dQk
Pk dPk
4 つの (1 次元) 調和振動子 :
=
,
= 2h3d Qk , k = 1, · · · , 4
ds
4
ds
1 P12 + P22 + P32 + P42
K2
エネルギーの値 : h3d =
− 2
8 Q21 + Q22 + Q23 + Q24
Q1 + Q22 + Q23 + Q24
ds
1
⇓ 時間変換 : s → t :
= 2
, k = 1, · · · , 4
dt
Q1 + Q22 + Q23 + Q24
KS 座標系上の 3D-Kepler :
dQk 1
Pk
dPk
2h2d Qk
=
,
=
, k = 1, · · · , 4,
dt
4 Q21 +Q22 +Q23 +Q24 dt
Q21 +Q22 +Q23 +Q24
エネルギーの値 : h3d
⇓ KS 変換 : x = Q21 −Q22 −Q23 +Q24 , y = 2(Q1 Q2 −Q3 Q4 ), z = 2(Q1 Q3 +Q2 Q4 ),
1 P1 Q1 −P2 Q2 −P3 Q3 +P4 Q4
1 P1 Q2 +P2 Q1 −P3 Q4 −P4 Q3
px =
,
p
=
,
y
2
Q21 +Q22 +Q23 +Q24
2
Q21 +Q22 +Q23 +Q24
1 P1 Q3 +P2 Q4 +P3 Q1 +P4 Q2
pz =
, P1 Q4 − P2 Q3 + P3 Q2 − P4 Q1 = 0
2
Q21 +Q22 +Q23 +Q24
13
3 次元 Kepler 運動と 4 重調和振動子 (2)
4 重調和振動子 + 1 拘束条件 ⇒ 3 次元 Kepler 問題 (3D-Kepler)
相対座標系上の 3D-Kepler :
dx
dpx
x
= px ,
= −K 2
3 ,
2
2
2
dt
dt
(x + y + z ) 2
dy
dpy
y
= py ,
= −K 2
3 ,
2
2
2
dt
dt
(x + y + z ) 2
dz
dpz
z
= pz ,
= −K 2
3 ,
2
2
2
dt
dt
(x + y + z ) 2
'
1& 2
K2
2
2
エネルギーの値 : h3d =
p + py + pz − !
2 x
x2 + y 2 + z 2
まとめ時間変換，正準変換によって，
“4 つの 1 次元の振動子 + 1 つの拘束条件からなる系” は 3D-Kepler
に変換される．
14
3 次元 Kepler 運動と 4 重調和振動子 (3)
1
保存量の関係
4 重調和振動子
相空間の次元 = 8 (自由度 4)
1
独立な保存量 (7 個)
1
2
調和振動子 1 ∼ 4 のエネルギー
(4 個)
Pk2
Hk :=
−h3d Q2k , k = 1,· · ·, 4
8
角運動量 (3 個)
J1 := P1 Q2 −P2 Q1
J2 := P1 Q4 −P4 Q1
J3 := P2 Q3 −P3 Q2 (= J2 と同じ値)
KS 座標系上の 3D-Kepler
相空間の次元 = 8 (自由度 4)
1
独立な保存量 (7 個)
1
エネルギー (1 個)
E=
2
H1 + H2 + H3 + H4 − 8K 2
Q21 + Q22
+ h3d = h3d
RL ベクトルの x, y 成分 (2 個)
ex := −8(H1 −H2 −H3 +H4), ey = F1 (H1 , ., H4 , J1 , ., J3)
3
角運動量 (3 個) hx = −J1 ,
hy = F2 (H1 , ., H4 , J1 , ., J3), hz = F3 (H1 , ., H4 , J1 , ., J3)
4
拘束条件 (1 個) 0 = J2 − J3
まとめ 4 重調和振動子がもつ独立の保存量 H1 , H2 , H3 , H4 , J1 , J2 , J3 が保
たれるならば，3D-Kepler がもつ全ての保存量も保たれる．
15
3 次元 Kepler 運動の全保存型差分 (1)
調和振動子 1,2,3,4 のエネルギー H1 , H2 ,H3 , H4 ; 時刻 s
⇓ 離散変分法 (3)
離散 4 重調和振動子
離散時刻 : s(j) , s(j+1)
(j+1)
(j)
(j+1)
(j)
(j+1)
(j)
(
)
Qk
−Qk
Pk
+Pk
Pk
−Pk
(j+1)
(j)
=
, (j+1)
= h2d Qk
+Qk , k = 1, .., 4; j = 0, 1, · · · ;
8
s(j+1) − s(j)
s
− s(j)
1
1
独立な保存量 : H1 , H2 , H3 , H4 , J1 , J2 , J3 の 7 個
⇓ 離散時間変換 s
1
(j)
離散 3D-Kepler
(j+1)
(j)
(j)
→t
s(j+1) −s(j)
1
: (j+1) (j) =
(j)
(j)
(j)
(j)
t
−t
(Q1 )2 +(Q2 )2 +(Q3 )2 +(Q4 )2
離散時刻 : t(n) , t(n+1)
(j+1)
(j)
(j+1)
(j)
(j+1)
(j)
Qk
−Qk
+Pk
Pk
−Pk
Qk
+Qk
1 Pk
=
,
=
h
, k = 1, .., 4; j = 0, 1, · · · ;
2d
(j+1)
(j)
(j)
(j+1)
(j)
(j)
2
2
8 |Q |
t
−t
t
−t
|Q |
(j)
(j)
(j)
(j)
|Q(j) |2 = (Q1 )2 +(Q2 )2 +(Q3 )2 +(Q4 )2 ;
1
独立な保存量 : E (エネルギー), ex , ey (RL ベクトルの 2 成分)，
hx , hy , hz (角運動量の 3 成分)，拘束条件の 7 個
16
3 次元 Kepler 運動の全保存型差分 (2)
離散 3D-Kepler において，
(運動の自由度)
1
=
=
(相空間の次元)
8
(独立な保存量の数)
7
−
−
が成り立つが，独立な保存量は全て，元の 3D-Kepler と同じもの．
=⇒ 離散 3D-Kepler の軌道は 3D-Kepler の軌道に一致する．
E=-0.5
E=-0.2
E= 0.0
E= 0.2
E= 0.5
z
0.5
0
-0.5
-1
-1.5
-2
-2.5
-3
0.5
-3.5
-4-3.5
0
-3
-0.5
-2.5
-2
-1.5
x
-1
-1
-0.5
0
y
-1.5
0.5
1 -2
3 次元離散 Kepler 問題による軌道の再現
17
離散 Kepler の時間補正 (1)
(仮定) LC (KS) 座標系で，
時刻 t0 における Kepler の位置 Q(t0 ) と
離散時刻 t(0) = t0 における離散 Kepler の位置 Q(0)
が一致している．
18
離散 Kepler の時間補正 (1)
(仮定) LC (KS) 座標系で，
時刻 t0 における Kepler の位置 Q(t0 ) と
離散時刻 t(0) = t0 における離散 Kepler の位置 Q(0)
が一致している．
19
離散 Kepler の時間補正 (2)
(仮定) LC (KS) 座標系で，
時刻 t0 における Kepler の位置 Q(t0 ) と
離散時刻 t(0) = t0 における離散 Kepler の位置 Q(0)
が一致している．
20
離散 Kepler の時間補正 (2)
(仮定) LC (KS) 座標系で，
時刻 t0 における Kepler の位置 Q(t0 ) と
離散時刻 t(0) = t0 における離散 Kepler の位置 Q(0)
が一致している．
21
離散 Kepler の時間補正 (3)
(仮定) LC 座標系上で，
t
(0)
t(1)
1
(0)
(0)
(0)
(0)
P 1 , P2 , Q 1 , Q 2
#
(1)
(1)
(1)
̸ t1 , P1 , P 2 , Q 1
=
&
= (P1 (t0 ), P2 (t0 ), Q1 (t0 ), Q2 (t0 ))
$ %
&
(1)
, Q2
= P1 (t1 ), P2 (t1 ), Q1 (t1 ) , Q2 (t1 )
時刻 t1 における 2D-Kepler 運動の座標
Q1 (t1 ) =
2
= t0 ,
%
'
t1
t0
dQ1 (t1 )
d
P1 (t1 )
Q1 (t) dt + Q1 (t0 ) ⇐⇒
=
dt
dt1
(Q1 (t1 )) 2 + (Q2 (t1 )) 2
離散時刻 t(1) における離散 2D-Kepler 運動の座標・運動量
⎧
⎪
P (1) = P (t0 ), Q(t0 ), t(1) の関数 ,
⎪
⎪
⎪
⎪
⎪
⎪
⎪
Q(1) = P (t0 ), Q(t0 ), t(1) の関数
⎪
⎪
⎨
離散 2D-Kepler ⇒
,
(1)
(0)
4
⎪
8h
(t
−
t
)|Q(t
)|
Q
(t
)
⎪
0
1 0
2d
⎪
%
&
⎪
⎪
(1)
4
(1)
(0)
2
2
⎪
dQ1
+2|Q(t0 )| 2|Q(t0 )| − h2d (t − t ) P1 (t0 )
⎪
⎪
⎪
⎪
=
.
/2
⎩ =⇒
dt(1)
2|Q(t0 )|4 − h2d (t(1) − t(0) )2
22
離散 Kepler の時間補正 (4)
1
時間補正の公式
離散時刻 t(1) で離散 2D-Kepler が与える位置 Q(1) に
2D-Kepler が到達する時刻 t1 が満たす関係式
(1)
(1)
dQ1
dt(1)
dQ1 (t1 ) dt1
=
⇐⇒
dt1 dt(1)
dQ1
dt1
dt(1)
=
dt(1)
dQ1 (t1 )
dt1
23
離散 Kepler の時間補正 (5)
1
離散 2D-Kepler における時間補正の公式 [楕円軌道 (h2d < 0) の場合]
4h2d (t(1) − t0 )|Q(t0 )|8 + 2(t(1) − t0 )|P (t0 )|2 |Q(t0 )|6 + 8|Q(t0 )|8 P (t0 ) · Q(t0 )
t1 = t0 +
.
/2
h2d 2|Q(t0 )|4 − h2d (t(1) − t0 )2
2h2d (t(1) − t0 )|Q(t0 )|4 + (t(1) − t0 )|P (t0 )|2 |Q(t0 )|2 + 8|Q(t0 )|4 P (t0 ) · Q(t0 )
.
/
−
h2d 2|Q(t0 )|4 − h2d (t(1) − t0 )2
!√
"
√
#
$
(1)
2
1
−2h2d (t − t0 )
+
(−h2d )−3/2
|P (t0 )|2 − h2d |Q(t0 )|2 tan−1
2
2
2|Q(t0 )|2
2
離散 3D-Kepler における時間補正の公式 [楕円軌道 (h3d < 0) の場合]
⎡
⎤
64h3d |Q(t0 )|6 + 8h23d (t(1) − t0 )2 |Q(t0 )|2
(t(1) − t0 )2 |Q(t0 )|2 ⎣ +8|P (t0 )|2 |Q(t0 )|4 + h3d (t(1) − t0 )2 |P (t0 )|2 ⎦
+32h3d (t(1) − t0 )|Q(t0 )|2 P (t0 ) · Q(t0 )
t1 = t0 +
.
/2
2h3d 8|Q(t0 )|4 − h3d (t(1) − t0 )2
!√
"
√ .
/
2
2
(1)
2 |P (t0 )| − 8h3d |Q(t0 )|
−2h3d (t − t0 )
−1
+
tan
4|Q(t0 )|2
8(−h3d )3/2
24
離散 Kepler の時間補正 (6)
2 次元離散 Kepler 問題の時間補正
25
離散 Kepler の時間補正 (7)
3 次元離散 Kepler 問題の時間補正
26
Results
Arbitrary time step size can be set in the discrete-time 2D and 3D-Kepler motions.
[ Property of discrete gradient method ]
1
1
2
Discrete-time 2D and 3D-Kepler motions preserve all conserved quantities .
[ Property of totally conservative integrator ]
−→ Discrete-time 2D and 3D-Kepler motions can exactly reproduce the orbits of the
original 2D and 3D-Kepler motions.
Discrete-time 2D and 3D-Kepler motions are explicit schemes .
−→ We can compute the time evolutions of the discrete-time 2D and 3D-Kepler motions
with high speed .
⇓
1
We presented time adjustment formulas which describe the time when the two and
three-dimensional Kepler motions arrive at the points located by the discrete Kepler
motions.
−→ Using a combination of discrete Kepler motions and the time adjustment formulas, we
can compute the time evolutions of 2D and 3D-Kepler motions precisely and with high
speed .
27

Download Report