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

RK4モジュールを使って落体の運動方程式。
せっかく計算本体を別モジュールにしたにも関わらず記述量が変わっていないのがなんとも。

program ode
  use runge_kutta_4
  implicit none
  double precision, parameter :: g = 9.80665      ! gravitational acceleration
  double precision, parameter :: t0 = 0           ! initial t
  double precision, parameter :: x0(2) = (/0, 0/) ! initial value at t0, 1st elem: position, 2nd elem: velocity
  double precision :: t
  double precision :: x(2)
  double precision, parameter :: udiv = 1         ! division unit
  integer, parameter :: ndiv = 1                  ! number of divisions per udiv
  double precision, parameter :: h = udiv / ndiv  ! delta-t
  integer, parameter :: niter = 10 / udiv * ndiv  ! number of iterations
  integer, parameter :: pdiv = ndiv               ! print interval
  integer :: n
  type (func_type), dimension(2) :: fns
  fns(1)%f => f1
  fns(2)%f => f2
  call set_func(fns)
  x = x0
  do n = 0, niter
    t = t0 + n * h
    if (mod(n, pdiv) == 0) print '(5F13.8)', t, x(1), exact1(t), x(2), exact2(t)
    x = next(t, x, h)
  end do
contains
  function f1(t, x)
    double precision :: f1
    double precision, intent(in) :: t
    double precision, dimension(:), intent(in) :: x
    f1 = x(2)
  end function f1
  function f2(t, x)
    double precision :: f2
    double precision, intent(in) :: t
    double precision, dimension(:), intent(in) :: x
    f2 = g
  end function f2
  double precision function exact1(t)
    double precision, intent(in) :: t
    exact1 = x0(1) + x0(2) * (t - t0) + g * (t - t0) ** 2 / 2
  end function exact1
  double precision function exact2(t)
    double precision, intent(in) :: t
    exact2 = x0(2) + g * (t - t0)
  end function exact2
end program ode

計算結果は当然のことながら前出のものに一致していたので省略。