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

落体の運動方程式をEuler法により数値的に解くFORTRANプログラム。

program ode
  implicit none
  double precision, parameter :: g = 9.80665      ! gravitational acceleration
  double precision, parameter :: t0 = 0           ! initial t
  double precision, parameter :: y0 = 0           ! initial position at t0
  double precision, parameter :: v0 = 0           ! initial velocity at t0
  double precision :: t
  double precision :: y = y0
  double precision :: v = v0
  double precision :: yn
  double precision, parameter :: udiv = 1         ! division unit
  integer, parameter :: ndiv = 1000               ! number of divisions per udiv
  double precision, parameter :: h = udiv / ndiv  ! delta-t
  integer, parameter :: niter = 2 / udiv * ndiv   ! number of iterations
  integer, parameter :: pdiv = ndiv / 10          ! print interval
  integer :: n
  do n = 0, niter
    t = t0 + n * h
    if (mod(n, pdiv) == 0) print '(5F13.8)', t, y, exacty(t), v, exactv(t)
    yn = nexty(t, y, v, h)
    v = nextv(t, y, v, h)
    y = yn
  end do
contains
  double precision function f1(t, y, v)
    double precision, intent(in) :: t, y, v
    f1 = v
  end function f1
  double precision function f2(t, y, v)
    double precision, intent(in) :: t, y, v
    f2 = g
  end function f2
  double precision function exacty(t)
    double precision, intent(in) :: t
    exacty = y0 + v0 * (t - t0) + g * (t - t0) ** 2 / 2
  end function exacty
  double precision function exactv(t)
    double precision, intent(in) :: t
    exactv = v0 + g * (t - t0)
  end function exactv
  double precision function nexty(t, y, v, h)
    double precision, intent(in) :: t, y, v, h
    double precision :: k1
    k1 = f1(t, y, v)
    nexty = y + h * k1
  end function nexty
  double precision function nextv(t, y, v, h)
    double precision, intent(in) :: t, y, v, h
    double precision :: k1
    k1 = f2(t, y, v)
    nextv = v + h * k1
  end function nextv
end program ode

の場合の結果は、

   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000
...snip
   1.00000000   4.89842176   4.90332508   9.80665016   9.80665016
...snip
   2.00000000  19.60349367  19.61330032  19.61330032  19.61330032

であり、落下距離の真値との誤差は、
初期時刻からの経過時間が1[s]で-0.00490[m]、2[s]で-0.00981[m]となり、
経過時間に比例していることが分かる。
また、刻み幅を半分にしたの場合、

   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000
...snip
   1.00000000   4.90087342   4.90332508   9.80665016   9.80665016
...snip
   2.00000000  19.60839700  19.61330032  19.61330032  19.61330032

誤差は経過時間1[s]で-0.00245[m]、2[s]で-0.00490[m]となり、
同一時刻で刻み幅に比例した大きさとなっている。