■人工衛星の軌道に関する運動方程式の数値解法

　人工衛星の軌道に関する運動方程式（２階の微分方程式）：
　　　x”＝　-GMx/r³
　　　y”＝　-GMy/r³
　　ここで、
　　　G　：　万有引力定数
　　　M　：　地球の質量
　　　r　：　地球の中心と人工衛星の距離、r² = x² + y²
を初期条件［時刻 t = 0 で x = r0, y = 0, 速度v = v0］のもとで解けば人工衛星の軌道が求められます。打上げ角度をαとすると、初期条件［時刻 t = 0 で速度v = v0］は［時刻 t = 0 で x' = v0sinα、y' = v0cosα］となります。


ここではこの微分方程式の近似解を代表的な数値解法であるオイラー（Euler）法およびルンゲ・クッタ（Runge-Kutta）法によって求める手順を示します。上記の微分方程式を次のように一般化して以後の説明をします。
　　x”＝　f (x, y)
　　y”＝　g (x, y)
　　初期条件：t=0 で x=x0, y=y0, x'=vx0, y'=vy0　　（注）x', y' は速度

それぞれの解法において、先ず1階の常微分方程式：x' = f (x,t) に対する手順を、次に上記の運動方程式に対する計算手順を記すことにします。

最も単純な解法ですが、一般に精度がよくありません。

[x' = f(x,t) に対する計算手順]
　時刻tにおける座標を(x)とし、時刻t+dtにおける座標を次の式で順次求めていきます。
・x  = x + f(x,t)*dt


[上記運動方程式に対する計算手順]
　時刻tにおける座標および速度を(x,y,vx,vy)とし、時刻t+dtにおける座標および速度を次の式で順次求めていきます。
・ax = f(x,y)
・vx = vx + ax*dt
・x  = x + vx*dt
・ay = g(x,y)
・vy = vy + ay*dt
・y  = y + vy*dt

計算精度がよく，また計算の安定性が高いのでよく使用されます。

[x' = f(x,t) に対する計算手順]
　時刻tにおける座標を(x)とし、時刻t+dtにおける座標を次の式で順次求めていきます。
・k1 = f(x, t)
・k2 = f(x+k1/2, t+dt/2)
・k3 = f(x+k2/2, t+dt/2)
・k4 = f(x+k3, t+dt)
・x = x + (k1 + 2*k2 * 2*k3 + k4)*dt/6
　すなわち、４個の値の加重平均を取ってｘの増分量を決定します。


[上記運動方程式に対する計算手順]
　先ず２階の連立微分方程式を次のように１階の連立微分方程式に変換します
（変位の微分が速度であることを利用します）。
　　x' ＝　vx
　　vx'＝　f(x, y)　　　（注）vx : ｘ方向速度
　　y' ＝　vy
　　vy'＝　g(x, y)　　　（注）vy : ｙ方向速度
すると、ルンゲ・クッタ法での計算手順は以下のようになります。
・kx1 = vx*dt
・kv1 = f(x, y)*dt
・ky1 = vy*dt
・kw1 = g(x, y)*dt

・kx2 = (vx+kv1/2)*dt
・kv2 = f(x+kx1/2, y+ky1/2)*dt
・ky2 = (vy+kw1/2)*dt
・kw2 = g(x+kx1/2, y+ky1/2)*dt

・kx3 = (vx+kv2/2)*dt
・kv3 = f(x+kx2/2, y+ky2/2)*dt
・ky3 = (vy+kw2/2)*dt
・kw3 = g(x+kx2/2, y+ky2/2)*dt

・kx4 = (vx+kv3)*dt
・kv4 = f(x+kx3, y+ky3)*dt
・ky4 = (vy+kw3)*dt
・kw4 = g(x+kx3, y+ky3)*dt

・x  = x  + (kx1 + 2*kx2 + 2*kx3 + kx4)/6
・vx = vx + (kv1 + 2*kv2 + 2*kv3 + kv4)/6
・y  = y  + (ky1 + 2*ky2 + 2*ky3 + ky4)/6
・vy = vy + (kw1 + 2*kw2 + 2*kw3 + kw4)/6

(注）前記運動方程式のルンゲ・クッタ法による計算手順として、下記のようにｘ成分計算後にｙ成分を
　　計算すると説明しているウェブサイトが見受けられますが、これでは精度がよくありません。
・kx1 = vx*dt
・kv1 = f(x, y)*dt
・kx2 = (vx+kv1/2)*dt
・kv2 = f(x+kx1/2, y)*dt
・kx3 = (vx+kv2/2)*dt
・kv3 = f(x+kx2/2, y)*dt
・kx4 = (vx+kv3)*dt
・kv4 = f(x+kx3, y)*dt
・x  = x  + (kx1 + 2*kx2 + 2*kx3 + kx4)/6
・vx = vx + (kv1 + 2*kv2 + 2*kv3 + kv4)/6
・同様にしてｙ成分についても計算。