\(\def\bm{\boldsymbol}\)\(\def\mb{\mathbf}\)\(\def\di{\displaystyle}\)\(\def\ve{\varepsilon_0}\)\(\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}\)\(\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}\)
1.制限３体問題
制限３体問題とは、2つの質点が円運動している重力場の中で、
第3の微小質点の運動を求めるものである。
2つの質点\(\,\mathrm{P}_1,\mathrm{P}_2\,\)の質量を\(\,M_1,M_2\,\)とし質点間の距離を\(\,a\,\)とすると、
ケプラーの第3法則から
\(\quad G(M_1+M_2)=\omega^2a^3\qquad\cdots\,\)(1)
が成り立つ。
1.1.回転座標系
\(\hspace{40mm}\)
\(\hspace{80mm}\)図1.回転座標系
この系の重心を座標原点にとり、回転していない座標系\(\,(O,XYZ)\,\)『慣性座標系』での第3の質点\(\,\mathrm{P}\,\)の座標を\(\,(X,Y,Z)\,\)として、回転座標系の変換を行う。
\(\,(X,Y)\,\)平面内で角速度\(\,\omega\,\)で回転している座標系を\(\,(O,xyz)\,\)とし、
その座標系での\(\,\mathrm{P}\,\)の座標を\(\,(x,y,z)\,\)とすると
\(\left. \begin{array}{l}
X=x\cos \omega t-y\sin \omega t \\[4pt]
Y=x\sin \omega t+y\cos \omega t
\end{array} \right\}\qquad\cdots\,\)(2)
である。ここで、(2)式を時間で微分する。
\(\quad\dot{X}=\dot{x}\cos\omega t-\omega x\sin\omega t-\dot{y}\sin\omega t-\omega y\cos\omega t\)
\(\hspace{9mm}=(\dot{x}-\omega y)\cos\omega t-(\dot{y}+\omega x)\sin\omega t\qquad\cdots\,\)(3)
\(\quad\ddot{X}=(\ddot{x}-\omega\dot{y})\cos\omega t-\omega(\dot{x}-\omega y)\sin\omega t\)
\(\hspace{25mm}-(\ddot{y}+\omega x)\sin\omega t-\omega(\dot{y}+\omega x)\cos\omega t\)
\(\hspace{9mm}=(\ddot{x}-2\omega\dot{y}-\omega^2x)\cos\omega t-(\ddot{y}+2\omega\dot{x}-\omega^2y)\sin\omega t\qquad\cdots\,\)(4)
\(\quad\dot{Y}=\dot{x}\sin\omega t+\omega x\cos\omega t+\dot{y}\cos\omega t-\omega y\sin\omega t\)
\(\hspace{9mm}=(\dot{x}-\omega y)\sin\omega t+(\dot{y}+\omega x)\cos\omega t\qquad\cdots\,\)(5)
\(\quad\ddot{Y}=(\ddot{x}-\omega\dot{y})\sin\omega t+\omega(\dot{x}-\omega y)\cos\omega t\)
\(\hspace{25mm}+(\ddot{y}-\omega\dot{x})\cos\omega t-\omega(\dot{y}+\omega x)\sin\omega t\)
\(\hspace{9mm}=(\ddot{x}-2\omega\dot{y}-\omega x)\sin\omega t+(\ddot{y}+2\omega\dot{x}-\omega^2y)\cos\omega t\qquad\cdots\,\)(6)
ここで、行列表示をすると(2)式は
\(\quad\left(\begin{matrix}X\\Y\end{matrix}\right)=\left(\begin{matrix}\cos\omega t&-\sin\omega t\\\sin\omega t&\cos\omega t\end{matrix}\right)\left(\begin{matrix}x\\y\end{matrix}\right)\qquad\cdots\,\)(7)
微分の式((3)式)と二次微分の式((4)式は、それぞれ次のように表す事が出来る
\(\quad\left(\begin{matrix}\dot{X}\\\dot{Y}\end{matrix}\right)=\left(\begin{matrix}\cos\omega t&-\sin\omega t\\\sin\omega t&\cos\omega t\end{matrix}\right)\left(\begin{matrix}\dot{x}-\omega y\\\dot{y}+\omega x\end{matrix}\right)\qquad\cdots\,\)(8)
\(\quad\left(\begin{matrix}\ddot{X}\\\ddot{Y}\end{matrix}\right)=\left(\begin{matrix}\cos\omega t&-\sin\omega t\\\sin\omega t&\cos\omega t\end{matrix}\right)\left(\begin{matrix}\ddot{x}-2\omega\dot{y}-\omega^2x\\\ddot{y}+2\omega\dot{x}-\omega^2y\end{matrix}\right)\qquad\dots\,\)(9)
1.2.計算モデル
\(\hspace{40mm}\)
慣性系の\(\,X\,\)軸上に、主星(質量\(=M_1\))と伴星(質量\(=M_2\))を置き、重心に原点を置く。
ここで、\(\dfrac{M_1}{M_1+M_2}=\mu_1\,,\quad \dfrac{M_2}{M_1+M_2}=\mu_2\,\,\)とする。
主星－伴星の距離を\(\,a\,\)とすると、\(M_1\,\)の座標は\(\,(-\mu_2a,0)\,\)、\(\,M_2\,\)の座標は\(\,(\mu_1a,0)\,\)となる。
点\(\,P\,\)の座標を\(\,(x,y)\,\)とすると\(\,\bm{r_1}\,,\,\bm{r_2}\,\)の長さは
\(\quad |\bm{r_1}|=\sqrt{\left\{-\mu_2a-x\right\}^2+y^2}\qquad\cdots\,\)(10)
\(\quad |\bm{r_2}|=\sqrt{\left\{\mu_1a-x\right\}^2+y^2}\qquad\cdots\,\)(11)
なお、\(P-G_c\,\)間の\(\,\bm{r}\,\)は\(\,\,|\bm{r}|=\sqrt{x^2+y^2}\,\,\)である。
1.3.運動方程式
点\(\,P\,\)の慣性系での座標を\(\,(X,Y)\,\)とし、質量を\(\,m\,\)とすると運動方程式は
\(\quad m\frac{d^2X}{dt^2}=m\ddot{X}=F_x\qquad\cdots\,\)(12)
\(\quad m\frac{d^2Y}{dt^2}=m\ddot{Y}=F_y\qquad\cdots\,\)(13)
となる。ここでベクトル量も(2)式に従い、座標変換すると
\(\quad \left. \begin{array}{c}
f_X=f_x\cos \omega t-f_y\sin nt \\[4pt]
f_Y=f_x\sin \omega t+f_y\cos \omega t\quad
\end{array} \right\} \qquad\cdots\,\)(15)
よって、回転座標系での運動方程式は
\(\quad \left. \begin{array}{l}
\ddot{x}-2\omega\dot{y}-\omega^2x=f_x\\[4pt]
\ddot{y}+2\omega\dot{x}-\omega^2y=f_y\quad\\[4pt]
\ddot{z}=f_z
\end{array} \right\} \qquad\cdots\,\)(17)
となる。ここで、点\(\,P\,\)に対する\(M_1\,,\,M_2\,\)からの万有引力は次のようになる。
\(\quad\dfrac{GM_1m\bm{r_1}}{r_1^3}+\dfrac{GM_2m\bm{r_2}}{r_2^3}\qquad\cdots\,\)(18)
1.4.無次元化
- 長さの単位として、両星間の距離\(\,\bm{a}\,\)を採用する。
- 質量の単位として、両星の質量の和に万有引力定数を掛けた\(\,\bm{G(M_1+M_2)}\,\)を採用する。
- 時間の単位として、以下の値を採用する
ケプラーの第３法則から\(\quad\dfrac{1}{\omega}=\sqrt{\dfrac{a^3}{G(M_1+M_2}}\,\,\)で次元は\(\,[T]\,\)である。
よって、\(\quad t=\sqrt{\dfrac{a^3}{G(M_1+M_2)}}\,\tau\quad\)とおけば\(\,\bm{\tau}\,\)が無次元の量になる。
(17)式と無次元化により運動方程式は以下のようになる。
\(\quad\left. \begin{array}{l}
\ddot{x}=2\omega\dot{y}+\omega^2x+f_x\\[4pt]
\ddot{y}=-2\omega\dot{x}+\omega^2y+f_y\quad\\[4pt]
\ddot{z}=f_z
\end{array} \right\}\qquad\cdots\,\)(20)
この微分は時間\(\,\bm{\tau}\,\)についてのものである。そして、\(\,f_x\,,\,f_y\,,\,f_z\,\)は、次のようになる。
\(\quad\left. \begin{array}{l}
f_x=-\dfrac{\mu_1(x-\mu_2)}{r_1^3}-\dfrac{\mu_2(x+\mu_1)}{r_2^3}\quad\\[4pt]
f_y=-\dfrac{\mu_1y}{r_1^3}-\dfrac{\mu_2 y}{r_2^3}\\[4pt]
f_z=-\dfrac{\mu_1z}{r_1^3}+\dfrac{\mu_2 z}{r_2^3}\\[4pt]
r_1^2=(x-\mu_2)^2+y^2+z^2\\[4pt]
r_2^2=(x+\mu_1)^2+y^2+z^2
\end{array} \right\}\qquad\cdots\,\)(22)
1.5.方程式の数値積分
方程式の積分で差分法によるものは、最初に差分表の作成が必要となる。それで、差分表を作らずに数値積分する方法として、ルンゲクッタ法がある。
ルンゲクッタ法はどのような常微分方程式にも使用できる。ただ、そのためには、方程式の系を１階の常微分方程式の系に直さなければならない。そのために、次のように連立させて解法する形とする。また変数は\(\,y_1\sim y_4\,\)のラベルとする。
\(\quad\left. \begin{array}{l}
y_1=x\,,\quad y_2=y\\[3pt]
y_3=\dot{x}\,,\quad y_4=\dot{y}\quad\\[3pt]
z=0\,,\quad \dot{z}=0
\end{array} \right\} \qquad\cdots\,\)(23)
新しい変数(\(\,y_1\sim y_4\,\))での方程式は次のようになる
\(\quad\left. \begin{array}{l}
\dfrac{dy_1}{dt}=y_3\\[6pt]
\dfrac{dy_2}{dt}=y_4\\[7pt]
\dfrac{dy_3}{dt}=2y_4+y_1+f_x\\[6pt]
\dfrac{dy_4}{dt}=-2y_3+y_2+f_y\\[6pt]
f_x=-\dfrac{\mu_1(y_1-\mu_2)}{r_1^3}-\dfrac{\mu_2(y_1+\mu_1)}{r_2^3}\quad\\[7pt]
f_y=-\dfrac{\mu_1y_2}{r_1^3}-\dfrac{\mu_2 y_2}{r_2^3}\\[7pt]
r_1^2=(y_1-\mu_2)^2+y_2^2\\[4pt]
r_2^2=(y_1+\mu_1)^2+y_2^2
\end{array} \right\} \)(25)
2.方程式の解法
初期値問題として数値的に解く。
\(\quad \dfrac{dy}{dx}=f(x,y)\qquad\cdots\,\)(30)
\(x=x_0\,\)の周りで Taylor 展開すると
\(\quad y(x_0+h)=y(x_0)+hy'(x_0)+\dfrac{h^2}{2!}y”(x_0)+\cdots\qquad\cdots\,\)(31)
ここで、\(y'(x_0)=f(x_0,y_0)\,\)なので
\(\quad y(x_0+h)=y(x_0)+hf(x_0,y_0)+O(h^2)\qquad\cdots\,\)(32)
となる。この解法には次に述べる Euler法などの手法がある
2.1 Euler 法
(32)式に従って計算するのが Euler法である。
Euler 法では、点\(\,(x_ i , y_ i)\,\)での微係数をそのまま\(\,x_{i+1}\,\)まで外挿して\(\,y_{i+1}\,\)を定める。
精度を上げるためさらに高次の Taylor展開を達成する。ここで次のような式で表すことを考える。
\(y(x_0+h)=y(x_0)+\mu_1hf(x_0,y_0)+\mu_2hf(x_0+\alpha h,y_0+\beta k)\qquad\cdots\,\)(33)
右辺第二項は、点\(\,(x_0,y_0)\,\)での傾き\(\,h・f(x_0,y_0)\,\)に\(\,\mu_1\,\)の重みを掛けたものであり、第三項は、\(\,x_0\,\)と\(\,x_0+h\,\)の間のある点\((x_0+\alpha h,y_0+\beta k)\,\)での傾き\(\,h・f(x_0+\alpha h,y_0+\beta k)\,\)に\(\,\mu_2\,\)の重みを掛けたものである。即ち、第二項と第三項とで、\(x_0\,\)と\(\,x_0+h\,\)の間の、中間の\(\,x\,\)における導関数を\(\,f(x,y)\,\)の関数値の適当な平均を用いて、\(y(x_0)\,\)から\(\,y(x_0+h)\,\)を求める事を考える。ここで、
\(\quad y(x_0+h)=y(x_0)+hf(x_0,y_0)+\dfrac{h^2}{2}\left(\pdr{f}{x}+f(x_0,y_0)\pdr{f}{y}\right)+\dfrac{h^3}{6}\left[\ppdr{f}{x^2}+\cdots\right.\qquad\cdots\,\)(34)
これを、\(xy\,\)2変数同時に展開して
\(\quad f(x_0+\alpha h,y_0+\beta k)=f_0+\alpha h\left(\pdr{f}{x}\right)_0+\beta k\left(\pdr{f}{y}\right)_0\qquad\cdots\,\)(35)
この(33)式と(35)式から、
\(\quad y(x_0+h)=y(x_0)+\mu_1 hf_0+\mu_2 h\left[f_0+\alpha\left(\pdr{f}{x}\right)_0+\beta\left(\pdr{f}{y}\right)_0\right]\\[-38pt]\)
\(\hspace{21mm}=y(x_0)+(\mu_1+\mu_2)hf_0+\alpha\mu_2h^2\left(\pdr{f}{x}\right)_0+\beta\mu_2hk\left(\pdr{f}{y}\right)_0\quad\cdots\,\)(36)
となる。この(36)式と(31)式を比べれば
\(\quad \mu_1+\mu_2=1\quad,\quad \alpha\mu_2=\beta\mu_2=\dfrac{1}{2}\qquad\cdots\,\)(37)
ここで、\(k=hf_0\,\)ならば、式(37)は、Taylor展開に\(\,O(h^2)\,\)の精度まで一致することが分かる。
特に、\(\alpha=\beta=1\quad,\quad\mu_1=\mu_2=\dfrac{1}{2}\,\)のときは
\(\quad y(x_0+h)+y(x_0)+\dfrac{1}{2}\left[hf(x_0,y_0)+hf(x_0+h,y_0+hf_0)\right]\qquad\cdots\,\)(38)
これは改良 Euler 公式と呼ばれる。
2.2 Runge-Kutta法
さらに高次の項まで展開すると以下のようになる。
\(\quad y(x_0+h)=y(x_0)+\dfrac{1}{6}\left\{k_1+2(k_2+k_3)+k_4\right\}\quad\cdots\,\)(40)
より以下の形になる。
\(\quad k_1=hf(x_0,y_0)\qquad\cdots\,\)(41)
\(\quad k_2=hf\left(x_0+\dfrac{h}{2},y_0+\dfrac{k_1}{2}\right)\qquad\cdots\,\)(42)
\(\quad k_3=hf\left(x_0+\dfrac{h}{2},y_0+\dfrac{k_2}{2}\right)\qquad\cdots\,\)(43)
\(\quad k_4=hf(x_0+h,y_0+k_3)\qquad\cdots\,\)(44)
(40)式は、Taylor 展開の四次の精度までを持つ。通常、この式を Runge-Kutta の公式と呼ぶ。
\(\hspace{20mm}\)
Runge-Kutta法では、まず点\(\,(x_i,y_i)\,\)での微係数から\(\,y\,\)の増分\(\,k_1\,\)を推定する。次に、点\(\,(x_i,y_i)\,\)といま得た点\(\,(x_{i+1},y_i+k1)\,\)の中間点での微係数から、これに平行な破線で示す直線によって増分\(\,k_2\,\)の推定を行う。これを更に繰り返して、第3の増分推定\(\,k_3\,\)を得る。
次いで、第3の推定点\(\,(x_i+\varDelta x,y_i+k_3)\,\)における微係数から第4の増分推定\(\,k_4\,\)を算出する。そして最後にこれら4つの増分推定を\(\,1:2:2:1\,\)の重みを付けて平均化して、これを最終的な増分推定として\(\,y_{i+1}\,\)を定めている。
2.3 Runge-Kutta-Gill法
(40)式を改良したのが Runge-Kutta-Gill法である。この方法は式(40)～(44)の係数を変化させたものである。
\(\quad y(x_0+h)=y(x_0)+\dfrac{1}{6}\left\{k_1+(2-\sqrt{2})k_2+(2+\sqrt{2})k_3+k_4\right\}\qquad\cdots\,\)(46)
\(\quad k_1=hf(x_0,y_0)\hspace{90mm}\cdots\,\)(47)
\(\quad k_2=hf\left(x_0+\dfrac{h}{2},y_0+\dfrac{k_1}{2}\right)\hspace{60mm}\cdots\,\)(48)
\(\quad k_3=hf\left[x_0+\dfrac{h}{2},y_0+(-1+\sqrt{2})+\left(1-\dfrac{1}{\sqrt{2}}\right)k_2\right]\qquad\cdots\,\)(49)
\(\quad k_4=hf(x_0+h,y_0+k_3)\hspace{50mm}\cdots\,\)(50)
2.4 Runge-Kutta-Gill法の実際の処理
手順は以下の通り
- 出発点の\(\,(x,y)\,\)を\(\,(x_0,y_0)\,\)とする\(\qquad\)初期値として\(\,q_0=0\,\)とする
- 定数として\(\,crk0=1-1/\sqrt{2}\quad,\quad crk1=1+1/\sqrt{2}\,\)とする
- \(y(x_0)=y_0\,\)をもとに\(\,k_1=h\cdot f(x_0,y_0)\,\)を計算\(\quad k_1=h\cdot f(x_0,y_0)\)
- \(k_1\,,\,y_0\,,q_0\,\)から\(\,y_1\,\)を計算する: 1回目\(\,[q_0=0]\,\)2回目\(\,[q_0=q_4]\)(前回の)\(\,y_1=y_0+(k_1-2q_0)/2\)
- \(k_1\,,\,q_0\,\)から\(\,q_1\,\)を計算する\(\quad q_1=q_0+3(k_1-2q_0)/2-k_1/2\)
- \(x_0+h/2\,\)と\(\,y1\,\)から\(\,k_2\,\)を計算する\(\quad k_2=h\cdot f(x_0+h/2,y_1)\)
- \(y_1\,,\,k_2\,,\,q_1\,\)から\(\,y_2\,\)を計算する\(\quad y_2=y_1+\,crk0\,(k_2-q_1)\)
- \(q_1\,,\,k_2\,\)から\(\,q_2\,\)を計算する\(\quad q_2=q_1+3\,crk0\,(k_2-q_1)-\,crk0\,k_2\)
- \(x_0+h/2\,\)と\(\,y2\,\)から\(\,k_3\,\)を計算する\(\quad k_3=h\cdot f(x_0+h/2,y_2)\)
- \(y_2\,,\,k_3\,,\,q_2\,\)から\(\,y_3\,\)を計算する\(\quad y_3=y_2+\,crk1\,(k_3-q_2)\)
- \(q_2\,,\,k_3\,q_3\,(\)前回の\(),\,k2\,\)から\(\,q_3\,\)を計算する\(\quad q_3=q_2+3\,crk1\,(k_3-q_3)-crk1\,k_2\)
- \(x_0+h\,\)と\(\,y3\,\)から\(\,k_4\,\)を計算する\(\quad k_4=h\cdot f(x_0+h,y_3)\)
- \(y_3\,,\,k_4\,,\,q_3\,\)から\(\,y_4\,\)を計算する\(\quad y_4=y_3+(k_4-2q_3)/6\)
- \(q_3\,,\,k_4\,\)から\(\,q_4\,\)を計算する\(\quad q_4=q_3+(k_4-2q_3)/2-k_4/2\quad\)次回の\(\,q_0\,\)になる
- (3)から(14)の処理を繰り返して、最終出力となる
3.ラグランジュ点
3体目を放出する点として、第1ラグランジュ点を採用する。
3.1 ラグランジュ点とは
ラグランジュ点とは, 質量差のある二つの天体が共通重心の周りをそれぞれ円軌道を描いて回っているとき, この二天体に比べて質量が無視できるほど小さな宇宙機をある速度を与えてこの軌道面内に置くと, 最初の二天体との相対位置を変えずに回り続けることができる位置のことであり, 五つの点が存在する。
3.2 第1ラグランジュ点 L\(_1\)
L\(_1\,\)は質量 M\(_1\,\)と M\(_2\,\)の互いに公転する2天体を結ぶ直線上で 2体の間に存在する。太陽系の場合、一般に地球よりも太陽に近い軌道を回る物体は地球よりも短い公転周期を持つ。
これは地球から重力で引かれる効果を無視した場合の話である。物体がちょうど地球と太陽の間にあると、地球の重力の効果によって、物体が太陽に引かれる力は弱められる。物体が地球に近ければ近いほどこの効果は大きい。この効果によって L\(_1\,\)では物体の公転周期が地球の公転周期と等しい。
3.3 L\(_1\,\)の計算
2天体間の距離を\(\,R_\,\)とする。天体の質量比 \(M_2/M_1\,\)を\(\,q\,\)とする。各々の天体から重心への距離をそれぞれ\(\,R_1\,,\,R_2\,\)とすると
\(\quad R_1=R\dfrac{M_2}{M_1+M_2}=R\dfrac{q}{1+q}\qquad\cdots\,\)(60)
\(\quad R_2=R\dfrac{M_1}{M_1+M_2}=R\dfrac{1}{1+q}\qquad\cdots\,\)(61)
\(\hspace{40mm}\)
\(M_1\,\)と\(\,M_2\,\)は互いの万有引力で引き合いながら共通重心\(\,G\,\)の周りを周回する。\(M_1\,\)の軌道が円軌道とすると、ケプラーの第3法則から
\(\quad M_1R_1\omega^2=M_1\dfrac{M_2}{M_1+M_2}R\omega^2=\dfrac{GM_1M_2}{R^2}\qquad\cdots\,\)(62)
L\(_1\,\)では両星からの万有引力と重心の周りを回転する遠心力が釣り合っている。L\(_1\,\)と重心\(\,G\,\)の距離を\(\,\ell_1\,\)L\(_1\,\)と\(\,M_2\,\)の距離を\(\,\ell_x\,\)とする。L\(_1\,\)点に質量\(\,m\,\)の物体があるとすると、それの運動方程式は
\(\quad\dfrac{GmM_1}{(R-\ell_1)^2}-\dfrac{GmM_2}{\ell_1^2}=m\omega^2(R_2-\ell_1)\)
\(\hspace{40mm}=mG\dfrac{M_1+M_2}{R^3}\left(R-\dfrac{M_2}{M_1+M_2}R-R\ell_x\right)\qquad\cdots\,\)(63)
ここで、式(63)の両辺を\(\,\dfrac{GmM_1}{R^2}\,\)で割り、\(M_2/M_1=q\,\)とおくと
\(\dfrac{R^2}{(R-\ell_1)^2}-\dfrac{R^2q}{\ell_1^2}=\dfrac{1+q}{R}\left(R-\dfrac{q}{1+q}R-\ell_1\right)\quad\cdots\,\)(65)
今回のプログラムでは、ここから5次方程式を立てて、Newton法で解くやり方と、両星の中間から引力と遠心力のバランスでの繰り返し計算の2通りで解いた。
4.運動軌跡
この連立方程式をルンゲクッタで解いたものが、下の図である。条件は伴星／主星の質量比を \(q=0.5\)とし、ラグランジュ点(L\(_1\,\))から物体を放した。物体が主星の周りでリサージュを描く様子が分かる。
\(\hspace{20mm}\)
座標系は主星伴星間の距離が1となっている。主星の半径が小さい(白色矮星)場合、このような形で降着円盤が形成されるが、アルゴル型などの場合、主星の半径が 0.1(比率) を超えるので、物体は回り込まずに中心からみて右上に衝突する形になる。
 |
 |
