万有引力

　2つの物体の間には,物体の質量に比例し,2物体間の距離の2乗に反比例する引力が作用する.

\[F = -G\frac{Mm}{r^2}\]

太陽の周りを回る地球の軌道を計算してみよう．

運動方程式

天体　\(i\)　が複数の天体　\(j\)　からの万有引力を受けて運動している場合，運動方程式は以下のように書ける．

(1)\[\frac{d^2 \boldsymbol{x}_i}{dt^2} = \sum_{j \neq i}Gm_j\frac{\boldsymbol{x}_j- \boldsymbol{x}_i}{|\boldsymbol{x}_j-\boldsymbol{x}_i|^3}\]

少々難しい書き方の気がするが，地球（\(i=e\))が太陽（\(j=s\)）の周りをまわっているとすれば，地球の運動方程式は，

\[\frac{d^2 \boldsymbol{x}_e}{dt^2} = Gm_s\frac{\boldsymbol{x}_s- \boldsymbol{x}_e}{|\boldsymbol{x}_s-\boldsymbol{x}_e|^3}\]

太陽を中心に固定し，　\(\boldsymbol{x}_e- \boldsymbol{x}_s = \boldsymbol{r}\)　と書けば，

\[\frac{d^2 \boldsymbol{r}}{dt^2} = - \frac{Gm_s}{|\boldsymbol{r}|^2} \frac{\boldsymbol{r}}{|\boldsymbol{r}|}\]

なじみの書き方になる． xy座標系で計算することにして，

\[\begin{split}\frac{d^2 x}{dt^2} = - \frac{Gm_s}{|\boldsymbol{r}|^2} \frac{x}{|\boldsymbol{r}|} \\ \frac{d^2 y}{dt^2} = - \frac{Gm_s}{|\boldsymbol{r}|^2} \frac{y}{|\boldsymbol{r}|}\end{split}\]

である．

無次元化（規格化？）

さて，運動方程式を数値計算すればよいのだが，万有引力定数 \(G = 6.67259 \times 10^{-11} \mathrm{m^3 s^{-2} kg^{-1}}\) ，太陽の質量 \(m_s = 1.989 \times 10^{30} \mathrm{kg}\) など非常に小さな数，大きな数が式中に現れている．これをそのまま計算するのは避けたい．

方程式を無次元化しよう．

まずは距離を \(x^* = x/a\)，時間を \(t^* = t/T\) のように無次元化する． \(a\) は，太陽と地球の距離（など）， \(T\)　は後で決める．

すると，

\[\frac{dx}{dt} = \frac{dx}{dt^*}\frac{dt^*}{dt}=\frac{1}{T}\frac{dx}{dt^*}=\frac{a}{T}\frac{dx^*}{dt^*}\]

同様に，

\[\frac{d^2x}{dt^2} =\frac{d}{dt}\left( \frac{a}{T}\frac{dx^*}{dt^*}\right)=\frac{d}{dt^*}\frac{dt^*}{dt}\left( \frac{a}{T}\frac{dx^*}{dt^*}\right)=\frac{a}{T^2}\frac{d^2 x^*}{d{t^*}^2}\]

これらを用いて運動方程式を書きなおすと，

\[\begin{split}\frac{a}{T^2}\frac{d^2x^*}{{dt^*}^2}=- \frac{Gm_s}{|a^2\boldsymbol{r}^*|^2} \frac{ax^*}{|a\boldsymbol{r}^*|} \\ \frac{a^3}{Gm_sT^2}\frac{d^2x^*}{{dt^*}^2}=- \frac{1}{|\boldsymbol{r}^*|^2} \frac{x^*}{|\boldsymbol{r}^*|}\end{split}\]

右辺からは単位をもった係数がなくなっている．距離がうまく無次元化できているならば大きさは1程度である．左辺の係数も大きさが1程度になるように選びたい．そこで，

\[\frac{a^3}{Gm_sT^2} = 1\]

になるように \(T\) を決める．すなわち，

\[T = \sqrt{\frac{a^3}{Gm_s}}\]

とする．すると結局運動方程式は，

\[\frac{d^2x^*}{{dt^*}^2}=- \frac{1}{|\boldsymbol{r}^*|^2} \frac{x^*}{|\boldsymbol{r}^*|}\]

となる．元の式と比べると，大きな（あるいは小さな）値をもった係数がなく簡潔である．

計算しよう

ボール投げの時と同様に，速度の変数を使う．（無次元変数を示す「*」は省略）

\[\begin{split}\dot x &= u, \dot u &= -\frac{x}{|\boldsymbol{r}^*|^3} \\ \dot y &= v, \dot v &= -\frac{y}{|\boldsymbol{r}^*|^3}\end{split}\]

ルンゲ－クッタ法を使って計算する．以下にコードを示す．

!+++++++++++++++++++++++++++++++++++++++++++++++
module globals

double precision,parameter::    G = 6.67384d-11
double precision                delta_t

end module globals
!+++++++++++++++++++++++++++++++++++++++++++++++
module subfunc

implicit none

contains
    function fx(x,y,u) result(xx)

        double precision,intent(in)::    x,y,u
        double precision        xx

        xx = u

    end function fx
    function fu(x,y,u) result(xx)

        double precision,intent(in)::    x,y,u
        double precision        r,xx

        r = sqrt(x**2+y**2)
        xx = -x / r**3

    end function fu
    function fy(x,y,u) result(xx)

        double precision,intent(in)::    x,y,u
        double precision        xx

        xx = u

    end function fy
    function fv(x,y,u) result(xx)

        double precision,intent(in)::    x,y,u
        double precision        r,xx

        r = sqrt(x**2+y**2)
        xx = -y / r**3

    end function fv
end module subfunc
!+++++++++++++++++++++++++++++++++++++++++++++++
program celestial

use globals
use subfunc

implicit none

integer                     i
integer,parameter::         n=100
double precision            tmax
double precision            t(n),x(n),y(n),u(n),v(n)
double precision            E(n)
double precision            kx1,kx2,kx3,kx4,ku1,ku2,ku3,ku4
double precision            ky1,ky2,ky3,ky4,kv1,kv2,kv3,kv4

x(1) = 1.0d0
y(1) = 0.0d0
u(1) = 0.0d0
v(1) = 1.0d0
t(1) = 0.0d0
tmax = 10.0d0

delta_t = (tmax-t(1))/dble(n)

do i = 1, n-1
    kx1 = fx(x(i),y(i),u(i))*delta_t
    ky1 = fy(x(i),y(i),v(i))*delta_t
    ku1 = fu(x(i),y(i),u(i))*delta_t
    kv1 = fv(x(i),y(i),v(i))*delta_t
    kx2 = fx(x(i)+0.5d0*kx1,y(i)+0.5d0*ky1,u(i)+0.5d0*ku1)*delta_t
    ky2 = fy(x(i)+0.5d0*kx1,y(i)+0.5d0*ky1,v(i)+0.5d0*kv1)*delta_t
    ku2 = fu(x(i)+0.5d0*kx1,y(i)+0.5d0*ky1,u(i)+0.5d0*ku1)*delta_t
    kv2 = fv(x(i)+0.5d0*kx1,y(i)+0.5d0*ky1,v(i)+0.5d0*kv1)*delta_t
    kx3 = fx(x(i)+0.5d0*kx2,y(i)+0.5d0*ky2,u(i)+0.5d0*ku2)*delta_t
    ky3 = fy(x(i)+0.5d0*kx2,y(i)+0.5d0*ky2,v(i)+0.5d0*kv2)*delta_t
    ku3 = fu(x(i)+0.5d0*kx2,y(i)+0.5d0*ky2,u(i)+0.5d0*ku2)*delta_t
    kv3 = fv(x(i)+0.5d0*kx2,y(i)+0.5d0*ky2,v(i)+0.5d0*kv2)*delta_t
    kx4 = fx(x(i)+kx3,y(i)+ky3,u(i)+ku3)*delta_t
    ky4 = fy(x(i)+kx3,y(i)+ky3,v(i)+kv3)*delta_t
    ku4 = fu(x(i)+kx3,y(i)+ky3,u(i)+ku3)*delta_t
    kv4 = fv(x(i)+kx3,y(i)+ky3,v(i)+kv3)*delta_t

    x(i+1) = x(i) + (kx1+2.0d0*kx2+2.0d0*kx3+kx4)/6.0d0
    u(i+1) = u(i) + (ku1+2.0d0*ku2+2.0d0*ku3+ku4)/6.0d0
    y(i+1) = y(i) + (ky1+2.0d0*ky2+2.0d0*ky3+ky4)/6.0d0
    v(i+1) = v(i) + (kv1+2.0d0*kv2+2.0d0*kv3+kv4)/6.0d0

    t(i+1) = t(i) + delta_t
end do

E = 0.5d0*(u**2 + v**2) - 1.0d0 / dsqrt(x**2+y**2)

open(unit=10,file='outputRK.dat',status='replace')
write(10,*)'t_rk, x_rk, y_rk, E_rk'
do i = 1,n
    write(10,*)t(i),',',x(i),',',y(i),',',E(i)
end do
close(10)

end program celestial
!+++++++++++++++++++++++++++++++++++++++++++++++

Warning

モジュールglobalはいらない．関数fuとfv内でのrの計算が重複している．関数の変数からtを（使わないので）削除した．またuだけでなくu,vとしたほうがよいのでは（使わないけど）．

座標(1.0,0.0)，速度(0.0,1.0)の初期条件で10秒(時間刻みは10/100=0.1)の計算を行っている． ボール投げ問題と比べて関数fxとかfuの中身が少々変わっているだけで，その他はほとんど変更なし．

軌道をプロットしてみると，円である．

この問題は解析的に軌道を求めることができる．極座標( \(r, \theta\) )表示で，初期条件が座標(1.0,0.0)，速度(0.0,v₀)の場合，

\[r = \frac{{v_0}^2}{1-(1-{v_0}^2)\cos\theta}\]

である．ちょうど \(|v_0|=1.0\) のとき円軌道， \(|v_0|=\sqrt{2}\) のとき放物線軌道， \(|v_0|>\sqrt{2}\) のとき双曲線軌道を描く．

せっかくなので，試してみよう．コード68行目「v(1) = 1.0d0」の値を変えればよい．

以下は， \(|v_0|= 0.8,1.0,1.2,\sqrt{2},1.6\) とした時の結果である．