Subsections

• 1 Euler 法
• 2 Runge-Kutta 法
• 3 例題
• 4 参考文献

1 Euler 法とRunge-Kutta 法

1 Euler 法

簡単のために、一元の微分方程式の場合を考えよう：

$${{dy}\over{dx}} = f(x,y) .$$ (20)

$x = x_0$の周りでの Taylor展開を利用する事はすぐに思いつくだろう：

$$y(x_0 + h) = y(x_0) + hy'(x_0) + {{h^2}\over{2!}}y''(x_0) + ... .$$ (21)

ここで、 $y'(x_0) = f(x_0, y_0)$ だから、$O(h)$までは、Taylor展開で求める事が 出来る。即ち、

$$y(x_0 + h) = y(x_0) + hf(x_0, y_0) + O(h^2)$$ (22)

誤差は

$${{h^2}\over{2!}}y''(x_0) = {{h^2}\over{2}} \left({{\partial f}\over{\partial x}} + f {{\partial f}\over{\partial y}}\right)$$ (23)

のオーダーである。この方法による計算は、Euler法と呼ばれる。

Figure 7.1: Euler法による積分の概念図。

2 Runge-Kutta 法

$x_0$ と $x_0+h$の間の、中間の$x$ における導関数、$f(x, y)$の関数値を何回か 計算し、それらの適当な平均を用いる事によって、高次の Taylor展開を達成しよう というのが、Runge-Kutta 法である。

次の例で説明しよう。

$$y(x_0+h) = y(x_0) + \mu_1 h f(x_0, y_0) + \mu_2 h f(x_0+\alpha h, y_0 + \beta k)$$ (24)

と表す。 右辺第二項は、点$(x_0, y_0)$での傾き $h\cdot f(x_0, y_0)$ に$\mu_1$ の 重みを掛けたものであり、第三項は、 $x_0$ と $x_0+h$の間のある点 $(x_0+\alpha h, y_0 + \beta k)$での傾き $h\cdot f(x_0+\alpha h, y_0 + \beta k)$ に$\mu_2$ の重みを掛けたものである。即ち、第二項と第三項とで、 $x_0$ と $x_0+h$の間の、中間の$x$ における導関数、$f(x, y)$の関数値 の適当な平均を用いて、$y(x_0)$ から$y(x_0+h)$を求める事を考える。 ここで、

$$f(x_0+\alpha h, y_0 + \beta k) = f_0 + \alpha h \left({{\p... ...t)_0 + \beta k \left({{\partial f}\over{\partial y}}\right)_0$$ (25)

だから、(7.6)式は、

$$y(x_0+h) = y(x_0) + \mu_1 h f_0 + \mu_2 h \left[ f_0 + \alpha h \left({{\par... ...l x}}\right)_0 + \beta k \left({{\partial f}\over{\partial y}}\right)_0 \right]$$ (26)

$$= y(x_0) + (\mu_1 + \mu_2) h f_0 + \alpha \mu_2 h^2 \left({{\partia... ...l x}}\right)_0 + \beta \mu_2 h k \left({{\partial f}\over{\partial y}}\right)_0$$

となる。 この式と(7.3)式とを比べれば、

$$\mu_1 + \mu_2 = 1$$ (27)

$$\alpha \mu_2 = \beta \mu_2 = 1/2$$ (28)

$$k = h f_0$$ (29)

ならば、(7.9)式は、Taylor展開に$O(h^2)$の精度まで一致する事が判る。 $\alpha$，$\beta$や重み$\mu_1$，$\mu_2$はユニークではない。

特に、 $\alpha = \beta =1$， $\mu_1 = \mu_2 =1/2$ の時は、

$$y(x_0+h) = y(x_0) + {{1}\over{2}}[h f_0 + h f_0(x_0 + h,y_0 + h f_0)]$$ (30)

で、これを改良Euler公式、或いは、Heun の二次公式と呼ぶ。 $\alpha = \beta =1/2$，$\mu_1=0$，$\mu_2=1$ の場合は、

$$y(x_0+h) = y(x_0) + h f(x_0 +h/2, y_0 + h/2 f_0 )$$ (31)

で、これは修正Euler公式と呼ばれる。 改良Euler公式も修正Euler公式も、幾何学的に考えれば、意味は明らかであろう。

Figure 7.2: 修正Euler法による積分の概念図。

ここでは証明は略すが、

$$y(x_0+h) = y(x_0) + {{1}\over{6}}\{k_1 + 2(k_2 + k_3) + k_4\}$$ (32)

但し、

$$k_1 = h f(x_0, y_0)$$ (33)

$$k_2 = h f(x_0 + h/2, y_0 + k_1/2)$$ (34)

$$k_3 = h f(x_0 + h/2, y_0 + k_2/2)$$ (35)

$$k_4 = h f(x_0 + h, y_0 + k_3)$$ (36)

とすると、 (7.14)式は、 Taylor展開の四次の精度までを持つ。通常、(7.14)式をRunge-Kuttaの公式と呼ぶ。

Runge-Kuttaの公式のサブルーチンプログラムの例を下に掲げる。このプログラムでは、 $n$元連立一次方程式で、 ${x_0, y_0, (dy/dx)_0, h}$は与えられているとして、$y(x_0+h)$ を計算する。 中間点での導関数を与えるプログラムは、別のサブルーチンで与えられているとして、 それをDERIVSと呼んでいる。この名前は仮の名前だから、引数として扱っている。 メインプログラムでは、external文を使って、導関数を与えるプログラムを指定 しておくだけで良い。

      subroutine RK4(y,dydx,n,x,h,yout,DERIVS)

      implicit NONE
      real*8 y,dydx,yout,yt,dyt,dym,
     &     x,xh,h,hh,h6
      integer nmax,n,i
      parameter (nmax=2)
      dimension y(n),dydx(n),yout(n),yt(nmax),dyt(nmax),dym(nmax)

      hh=h*0.5
      h6=h/6.
      xh=x+hh
      Do i=1,n
        yt(i)=y(i)+hh*dydx(i)
      End do
      call DERIVS(xh,yt,dyt)
      Do i=1,n
        yt(i)=y(i)+hh*dyt(i)
      End do
      call DERIVS(xh,yt,dym)
      Do i=1,n
        yt(i)=y(i)+h*dym(i)
        dym(i)=dyt(i)+dym(i)
      End do
      call DERIVS(x+h,yt,dyt)
      Do i=1,n
        yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
      End do

      Return
      End

3 例題

Runge-Kutta法とEuler法と修正Euler法を使って、次の連立一階常微分方程式を解き、 誤差を評価せよ。

$${{d^2 y}\over{dx^2}} + y = 0 .$$ (37)

但し、初期値は $x=0$で$y=1$，$dy/dx=0$ とし、 区間 $x= [0, \pi/2]$ について解け。 尚、先に挙げたサブルーチン RK4 の使用例として、 Runge-Kutta法でのサンプルプログラムを下に掲げる。

      program main

      implicit NONE
      real*8 y,dydx,x,h,yout,pi
      integer i
      parameter(i=2,h=0.001,pi=3.141593)
      dimension y(i),dydx(i),yout(i)
      external diff_eq

      open(unit=1,file='RK_sample1.dat',status='new')
c initial value
      x    = 0.
      y(1) = 1.
      y(2) = 0. 
c do loop
      do while (x.le.pi/2.)
        write(1,100) x,y(1),y(2)       
        call diff_eq(x,y,dydx)
        call rk4(y,dydx,2,x,h,yout,diff_eq)
        y(1) = yout(1)
        y(2) = yout(2)
        x    = x + h
      end do
100   format(11x,f6.3,1p2e12.4)
      close(1)
      stop
      end


      subroutine diff_eq(x,y,dydx)

      implicit NONE
      real*8 y,dydx,x,n
      dimension y(2),dydx(2)

      dydx(1)=y(2)
      dydx(2)=-y(1)
      return
      end

4 参考文献

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. 1992, Numerical Recipes in FORTRAN 77 (Second Edition), (Cambridge University Press, Cambridge).

http://cfatab.harvard.edu/nr/nronline.html