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

Heun法を使ってみる。

...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, x + h * k1)
    nextx = x + h / 2 * (k1 + k2)
  end function nextx
...snip

計算してみると、

   0.00000000   0.00000000   0.00000000
   0.10000000   0.10033469   0.10033467
   0.20000000   0.20271007   0.20271004
   0.30000000   0.30933630   0.30933625
   0.40000000   0.42279328   0.42279322
   0.50000000   0.54630257   0.54630249
   0.60000000   0.68413689   0.68413681
   0.70000000   0.84228846   0.84228838
   0.80000000   1.02963859   1.02963856
   0.90000000   1.26015814   1.26015822
   1.00000000   1.55740734   1.55740772
   1.10000000   1.96475845   1.96475966
   1.20000000   2.57214793   2.57215162
   1.30000000   3.60208942   3.60210245
   1.40000000   5.79781591   5.79788372
   1.50000000  14.10022034  14.10141995

t=1.5で9e-5弱の相対誤差であり中点法と同程度の誤差を持つ。