微分方程式の数値解析 #9

FORTRAN

太陽を周回する地球の運動について見てみる。
解くべき運動方程式は時刻tでの軌道面における地球の直交座標を(x(t),y(t))として、
\begin{eqnarray}\ddot{x}=-GMx(x^2+y^2)^{-3/2}\\ \ddot{y}=-GMy(x^2+y^2)^{-3/2}\end{eqnarray}
である。太陽は原点から動かないものとする。GMは万有引力定数と太陽質量の積である。

program ode
  use runge_kutta_4
  implicit none
  double precision, parameter :: gm = 9.9069307177135d8 ! prod of gravitational constant and sun mass [(10^4km)^3/day^2]
  double precision, parameter :: t0 = 0                 ! initial t
  double precision, parameter :: x0(4) = (/1.49597871d4, 0d0, 0d0, 2.572992d2/) ! (x0, y0, vx0, vy0)
  double precision :: t
  double precision :: x(4)
  double precision, parameter :: udiv = 1         ! division unit
  integer, parameter :: ndiv = 4                  ! number of divisions per udiv
  double precision, parameter :: h = udiv / ndiv  ! delta-t
  integer, parameter :: niter = 366 / udiv * ndiv ! number of iterations
  integer, parameter :: pdiv = ndiv               ! print interval
  integer :: n
  type (func_type), dimension(4) :: fns
  fns(1)%f => f1
  fns(2)%f => f2
  fns(3)%f => f3
  fns(4)%f => f4
  call set_func(fns)
  x = x0
  do n = 0, niter
    t = t0 + n * h
    if (mod(n, pdiv) == 0) print '(5F13.5)', t, x(1), x(2), x(3), x(4)
    x = next(t, x, h)
  end do
contains
  double precision function f1(t, x)
    double precision, intent(in) :: t
    double precision, dimension(:), intent(in) :: x
    f1 = x(3)
  end function f1
  double precision function f2(t, x)
    double precision, intent(in) :: t
    double precision, dimension(:), intent(in) :: x
    f2 = x(4)
  end function f2
  double precision function f3(t, x)
    double precision, intent(in) :: t
    double precision, dimension(:), intent(in) :: x
    f3 = -gm * x(1) * (X(1) ** 2 + X(2) ** 2) ** (-1.5)
  end function f3
  double precision function f4(t, x)
    double precision, intent(in) :: t
    double precision, dimension(:), intent(in) :: x
    f4 =  -gm * x(2) * (X(1) ** 2 + X(2) ** 2) ** (-1.5)
  end function f4
end program ode

結果は、経過時間/day、軌道面でのX座標/10^4km、Y座標/10^4km、X方向の速度/(10^4km/day)、Y方向の速度/(10^4km/day)で、

      0.00000  14959.78710      0.00000      0.00000    257.29920
      1.00000  14957.57376    257.28651     -4.42657    257.26113
      2.00000  14950.93441    514.49689     -8.85182    257.14694
...snip
    364.00000  14957.18459   -278.98864      4.79995    257.25444
    365.00000  14959.77135    -21.70561      0.37344    257.29893
    366.00000  14957.93144    235.58384     -4.05318    257.26728

であり、刻み幅0.25dayとあまり細かくないにも関わらず365日強で公転する様が再現される。
刻み幅1dayでも以下のようになる。

...snip
    364.00000  14957.18459   -278.98843      4.79994    257.25444
    365.00000  14959.77134    -21.70540      0.37344    257.29893
    366.00000  14957.93143    235.58405     -4.05318    257.26728

刻み幅を1e-3dayにした結果を以下に示す。

...snip
    364.00000  14957.18459   -278.98864      4.79995    257.25444
    365.00000  14959.77135    -21.70561      0.37344    257.29893
    366.00000  14957.93144    235.58384     -4.05318    257.26728

ここまで細かくしなくてもいいことが分かる。