微分方程式の数値解析 #1 - 借りてきた猫のように静か

関数 $x=x(t)$ についての常微分方程式の初期値問題
$\begin{eqnarray}x'&=&x^2+1\\x(0)&=&0\end{eqnarray}$
の解析解は $x(t)=\tan t$ である。
これをEuler法で数値解析するプログラムをFORTRAN、古流書式でなく90/95で書いてみる。
微分方程式の初期値問題が
$\begin{eqnarray}x'&=&f(t,x)\\x(t_0)&=&x_0\end{eqnarray}$
の形の場合、刻み幅を $h$ とし、
$\begin{eqnarray}t_n&=&t_0+nh\\x_n&=&x(t_n)\end{eqnarray}$
と定義すると、Euler法は、
$\begin{eqnarray}x_{n+1}&=&g(t_n,x_n,h)\\g(t_n,x_n,h)&=&x_n+hf(t_n,x_n)\end{eqnarray}$
と書けるので、素直に実装する。

! Euler's method
program ode
  implicit none
  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 = 1.5 / udiv * ndiv ! number of iterations
  double precision, parameter :: t0 = 0           ! initial t
  double precision, parameter :: x0 = 0           ! initial value at t0
  integer, parameter :: pdiv = ndiv / 10          ! print interval
  double precision :: x = x0
  double precision :: t
  integer :: n
  do n = 0, niter
    t = t0 + n * h
    if (mod(n, pdiv) == 0) print '(3F13.8)', t, x, exactx(t)
    x = nextx(t, x, h)
  end do
contains
  double precision function f(t, x)
    double precision, intent(in) :: t, x
    f = x**2 + 1
  end function f
  double precision function exactx(t)
    double precision, intent(in) :: t
    exactx = tan(t)
  end function exactx
  double precision function nextx(t, x, h)
    double precision, intent(in) :: t, x, h
    nextx = x + h * f(t, x)
  end function nextx
end program ode

実行すると、

   0.00000000   0.00000000   0.00000000
   0.10000000   0.10032963   0.10033467
   0.20000000   0.20268911   0.20271004
   0.30000000   0.30928626   0.30933625
   0.40000000   0.42269643   0.42279322
   0.50000000   0.54613318   0.54630249
   0.60000000   0.68385543   0.68413681
   0.70000000   0.84183091   0.84228838
   0.80000000   1.02889559   1.02963856
   0.90000000   1.25893091   1.26015822
   1.00000000   1.55530569   1.55740772
   1.10000000   1.96093380   1.96475966
   1.20000000   2.56446697   2.57215162
   1.30000000   3.58384542   3.60210245
   1.40000000   5.73758911   5.79788372
   1.50000000  13.59898621  14.10141995

となり、刻み幅1e-3の場合、t=1.5で厳密解との相対誤差4e-2弱で求められた。