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

Euler法よりもマシなはずの中点法で計算してみる。
Euler法との違いは漸化式のの定義部分だけである。

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

実行結果は、

   0.00000000   0.00000000   0.00000000
   0.10000000   0.10033466   0.10033467
   0.20000000   0.20271002   0.20271004
   0.30000000   0.30933621   0.30933625
   0.40000000   0.42279316   0.42279322
   0.50000000   0.54630239   0.54630249
   0.60000000   0.68413664   0.68413681
   0.70000000   0.84228810   0.84228838
   0.80000000   1.02963807   1.02963856
   0.90000000   1.26015733   1.26015822
   1.00000000   1.55740601   1.55740772
   1.10000000   1.96475607   1.96475966
   1.20000000   2.57214305   2.57215162
   1.30000000   3.60207690   3.60210245
   1.40000000   5.79776617   5.79788372
   1.50000000  14.09953074  14.10141995

Euler法と同じ刻み幅で、t=1.5での相対誤差が1e-4強でありかなり改善されている。