C.3 Runge-Kutta (ルンゲ-クッタ)法

前節で解説した Euler 法は簡単で、これですべてが片付けば喜ばしいので すが、残念ながらあまり効率が良くありません。高精度の解を計算するために は、分割数 $N$ をかなり大きく取る(=大量の計算をする)必要があります。特 別な場合^C.9を除けば、実際に使われることは滅多にないでしょう。率直 にいって実用性は低いです。

より高精度の公式は、現在まで様々なものが開発されていますが、比較的簡 単で、精度がまあまあ高いものに Runge-Kutta 法と呼ばれるものがあります。 それは $x_j$ から $x_{j+1}$ を求める漸化式として次のものを用います。

$$k_1 = h f(t_j,x_j)$$

$$k_2 = h f(t_j+h/2,x_j+k_1/2)$$

$$k_3 = h f(t_j+h/2,x_j+k_2/2)$$

$$k_4 = h f(t_j+h,x_j+k_3)$$

$$x_{j+1} = x_j+\frac{1}{6}(k_1+2k_2+2k_3+k_4)$$

これがどうやって導かれたものかは解説しません。まずは使ってみましょう。

* reidai5-3.f -- 微分方程式の初期値問題を Runge-Kutta 法で解く
      program rei53
*       開始時刻と終了時刻
      real a,b
      parameter (a=0.0,b=1.0)
*       変数と関数の宣言
      integer N,j
      real t,x,h,x0,f,k1,k2,k3,k4
      external f
*       初期値
      x0 = 1.0
*      区間の分割数 N を入力してもらう
      write(*,*) ' N='
      read(*,*) N
*      小区間の幅
      h=(b-a)/N
*      開始時刻と初期値のセット
      t=a
      x=x0
      write(*,*) t,x
*      Runge-Kutta 法による計算
      do j=0,N-1
         k1 = h * f(t, x)
         k2 = h * f(t + h / 2, x + k1 / 2)
         k3 = h * f(t + h / 2, x + k2 / 2)
         k4 = h * f(t + h, x + k3)
         x = x + (k1 + 2 * k2 + 2 * k3 + k4) / 6
         t = t + h
         write(*,*) t,x
      end do
      end
**********************************************************************
*       微分方程式 x'=f(t,x) の右辺の関数 f の定義 --- これも前と同じだから略

コンパイル・実行の結果は次のようになるはずです。

waltz11% f77o reidai5-3.f
Fortran 90 diagnostic messages: program name(f)
  jwd2008i-i  ``reidai5-3.f'', line 30: この仮引数は,副プログラム中で使用されていま
せん.(名前:t)
waltz11% reidai5-3
  N=
10
 0.000000000e+00 1.00000000
 0.100000001 1.10517085
 0.200000003 1.22140265
 0.300000012 1.34985852
 0.400000006 1.49182427
 0.500000000 1.64872062
 0.600000024 1.82211792
 0.700000048 2.01375151
 0.800000072 2.22553945
 0.900000095 2.45960140
 1.00000012 2.71827984
waltz11%

たった $10$ 等分なのに相対誤差が $10^{-6}$ 以下になっています ($\because$ $x(1)=e^1=e=2.7182818284\cdots$ であるので、相対誤差 = $\vert e-2.71827984\vert / e=7.31\cdots\times 10^{-7}$)。 Runge-Kutta 法の公 式は Euler 法よりは大部面倒ですが、それに見合うだけの価値があることが 納得できるでしょう。