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

定番のclassical 4th-order Runge-Kutta法、略してRK4を使ってみる。

...snip
  double precision function nextx(t, x, h)
    double precision, intent(in) :: t, x, h
    double precision :: k1, k2, k3, k4
    k1 = f(t, x)
    k2 = f(t + h / 2, x + h / 2 * k1)
    k3 = f(t + h / 2, x + h / 2 * k2)
    k4 = f(t + h, x + h * k3)
    nextx = x + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
  end function nextx
...snip

結果は、

   0.00000000   0.00000000   0.00000000
   0.10000000   0.10033467   0.10033467
   0.20000000   0.20271004   0.20271004
   0.30000000   0.30933625   0.30933625
   0.40000000   0.42279322   0.42279322
   0.50000000   0.54630249   0.54630249
   0.60000000   0.68413681   0.68413681
   0.70000000   0.84228838   0.84228838
   0.80000000   1.02963856   1.02963856
   0.90000000   1.26015822   1.26015822
   1.00000000   1.55740772   1.55740772
   1.10000000   1.96475966   1.96475966
   1.20000000   2.57215162   2.57215162
   1.30000000   3.60210245   3.60210245
   1.40000000   5.79788372   5.79788372
   1.50000000  14.10141994  14.10141995

t=1.5での厳密解との相対誤差は1e-9以下である。
刻み幅はそれほど細かくないにも関わらず良い結果が出ている。