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

FORTRAN

落体の運動方程式をRK4により数値的に解くFORTRANプログラム。
Euler法によるプログラムのうちy(t)とv(t)からy(t+h)とv(t+h)を求める漸化式部分をRK4に変更しただけ。
同じ計算をnextyとnextvで行っているのでおいしくないプログラムになっている。
もっとまとめられて綺麗になるはずだが面倒なのでこのままで。

...snip
  double precision function nexty(t, y, v, h)
    double precision, intent(in) :: t, y, v, h
    double precision :: ky1, ky2, ky3, ky4, kv1, kv2, kv3, kv4
    ky1 = f1(t, y, v)
    kv1 = f2(t, y, v)
    ky2 = f1(t + h / 2, y + h / 2 * ky1, v + h / 2 * kv1)
    kv2 = f2(t + h / 2, y + h / 2 * ky1, v + h / 2 * kv1)
    ky3 = f1(t + h / 2, y + h / 2 * ky2, v + h / 2 * kv2)
    kv3 = f2(t + h / 2, y + h / 2 * ky2, v + h / 2 * kv2)
    ky4 = f1(t + h, y + h * ky3, v + h * kv3)
    kv4 = f2(t + h, y + h * ky3, v + h * kv3)
    nexty = y + h / 6 * (ky1 + 2 * ky2 + 2 * ky3 + ky4)
  end function nexty
  double precision function nextv(t, y, v, h)
    double precision, intent(in) :: t, y, v, h
    double precision :: ky1, ky2, ky3, ky4, kv1, kv2, kv3, kv4
    ky1 = f1(t, y, v)
    kv1 = f2(t, y, v)
    ky2 = f1(t + h / 2, y + h / 2 * ky1, v + h / 2 * kv1)
    kv2 = f2(t + h / 2, y + h / 2 * ky1, v + h / 2 * kv1)
    ky3 = f1(t + h / 2, y + h / 2 * ky2, v + h / 2 * kv2)
    kv3 = f2(t + h / 2, y + h / 2 * ky2, v + h / 2 * kv2)
    ky4 = f1(t + h, y + h * ky3, v + h * kv3)
    kv4 = f2(t + h, y + h * ky3, v + h * kv3)
    nextv = v + h / 6 * (kv1 + 2 * kv2 + 2 * kv3 + kv4)
  end function nextv
...snip

t_0=0,y_0=0,v_0=0,h=1e-3での結果は、

   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000
...snip
   1.00000000   4.90332508   4.90332508   9.80665016   9.80665016
...snip
   2.00000000  19.61330032  19.61330032  19.61330032  19.61330032

であり、落下距離も落下速度も真値に一致する。

ところで、この方程式についてはRK4では時間の刻み幅に無関係に誤差は無いわけであるから、
繰り返し計算による計算誤差の累積を考えれば、
刻み幅は必要な時刻の値だけが得られる程度に大きい方がよいことになる。
t_0=0,y_0=0,v_0=0,h=1で10[s]まで計算すると、

   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000
   1.00000000   4.90332508   4.90332508   9.80665016   9.80665016
   2.00000000  19.61330032  19.61330032  19.61330032  19.61330032
   3.00000000  44.12992573  44.12992573  29.41995049  29.41995049
   4.00000000  78.45320129  78.45320129  39.22660065  39.22660065
   5.00000000 122.58312702 122.58312702  49.03325081  49.03325081
   6.00000000 176.51970291 176.51970291  58.83990097  58.83990097
   7.00000000 240.26292896 240.26292896  68.64655113  68.64655113
   8.00000000 313.81280518 313.81280518  78.45320129  78.45320129
   9.00000000 397.16933155 397.16933155  88.25985146  88.25985146
  10.00000000 490.33250809 490.33250809  98.06650162  98.06650162

となる。
とはいえ、刻み幅1e-3[s]程度では10[s]の値はこの1[s]の時と変わらず誤差は見られなかった。
刻み幅1e-6[s]でやっと

...snip
  10.00000000 490.33250806 490.33250809  98.06650160  98.06650162

この程度の誤差が現れる。
累積誤差だけでなく計算時間的な問題もあるので、
刻み幅が大きくても構わないのであれば大きくすべきであろう。