微分方程式の数値解析 #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
計算結果は当然のことながら前出のものに一致していたので省略。