4次ルンゲ・クッタ法

微分方程式なら任せろーバリバリバリー

ルンゲクッタ法ってもともとどういうもの?理論は?刻み幅\(h\)を自動的に制御する方法について知りたい方は、ルンゲクッタ法の系統的扱いと刻み幅制御へどうぞ。

4次ルンゲ・クッタ法は微分方程式を数値的に解く手段です。

ルンゲクッタ法が良く使われる理由は、ひとえにプログラムの実装のしやすさです。ルンゲクッタ法は非常に良い!というアルゴリズムではありませんが、他の方法よりもシンプルで、プログラムに組み込みやすいのです。

”4次”は問題の解の関数をテイラー展開した場合、4次までは一致するように作られた、という意味です。
例えば微分方程式
\(
\displaystyle \frac{dy}{dx}=g(x,y)
\)

を考えます。
数値計算では初期値(例えばx=0の時、y=1など)を決めて、そこからxをhだけ増やし、微分方程式というルールに従って関数\(y(0+h),y(0+h+h),y(0+h+h+h),\cdots\)を作り上げていきます。
この時、4次ルンゲクッタ法で求められる答えの関数というのは
\(
\displaystyle y(x+h)=y(x)+hg(x,y)+\frac{h^2}{2!}\frac{dg}{dx}+\frac{h^3}{3!}\frac{d^2g}{dx^2}+\frac{h^4}{4!}\frac{d^3g}{dx^3}+O(h^5)
\)

という関数になります。
\(
g(x,y)=-xy
\)

(ただし、初期条件\(x=0\)で\(y=1\))
である場合、微分方程式の解析解は
\(
\displaystyle y=e^{-x^2/2}
\)

であるため、4次ルンゲクッタ法によって導かれる答えは、
\(
\displaystyle y(x+h)=y(x)\left[1-hx+\frac{h^2}{2}(x^2-1)+\frac{h^3}{6}(-x^3+3x)+\frac{h^4}{24}(x^4-6x^2+3)\right]+O(h^5)
\)

となります。

本題に入りましょう。
4次ルンゲクッタ法は6つのステップが必要となります。
初期値を\((x_0,y_0)\)と書くと、(\(y_0=y(x_0)\)です。)

  1. \((x_0,y_0)\)より\(k_a\)を求める
  2. \(k_a\)より\(k_b\)を求める
  3. \(k_b\)より\(k_c\)を求める
  4. \(k_c\)より\(k_d\)を求める
  5. \(k_a,k_b,k_c,k_d\)より\(y(x_0+h)\)を求める
  6. \((x_0+h,y(x_0+h))\)を初期値だと思って手順1に戻る。

という感じです。
3章:連立ルンゲ・クッタ法による微分方程式の解を参考にすると、数値計算での4次ルンゲクッタ法の計算スキームは以下のようになります。

解きたい微分方程式を連立1次微分方程式の形で書くと一般的にはこう書けます。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{dy_1}{dx} &= f_1(x,y_1,y_2,\cdots,y_N) \\
\frac{dy_2}{dx} &= f_2(x,y_1,y_2,\cdots,y_N) \\
&\vdots \\
\frac{dy_N}{dx} &= f_N(x,y_1,y_2,\cdots,y_N) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
まず係数\(k_a\)を求めます。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{a1} &= hf_1(x,y_1,y_2,\cdots,y_N) \\
k_{a2} &= hf_2(x,y_1,y_2,\cdots,y_N) \\
&\vdots \\
k_{aN} &= hf_N(x,y_1,y_2,\cdots,y_N) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
次に\(k_b\)。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{b1} &= hf_1(x+\frac{h}{2},y_1+\frac{k_{a1}}{2},y_2+\frac{k_{a2}}{2},\cdots,y_N+\frac{k_{aN}}{2}) \\
k_{b2} &= hf_2(x+\frac{h}{2},y_1+\frac{k_{a1}}{2},y_2+\frac{k_{a2}}{2},\cdots,y_N+\frac{k_{aN}}{2}) \\
&\vdots \\
k_{bN} &= hf_N(x+\frac{h}{2},y_1+\frac{k_{a1}}{2},y_2+\frac{k_{a2}}{2},\cdots,y_N+\frac{k_{aN}}{2}) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
そして\(k_c\)。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{c1} &= hf_1(x+\frac{h}{2},y_1+\frac{k_{b1}}{2},y_2+\frac{k_{b2}}{2},\cdots,y_N+\frac{k_{bN}}{2}) \\
k_{c2} &= hf_2(x+\frac{h}{2},y_1+\frac{k_{b1}}{2},y_2+\frac{k_{b2}}{2},\cdots,y_N+\frac{k_{bN}}{2}) \\
&\vdots \\
k_{cN} &= hf_N(x+\frac{h}{2},y_1+\frac{k_{b1}}{2},y_2+\frac{k_{b2}}{2},\cdots,y_N+\frac{k_{bN}}{2}) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
最後に\(k_d\)。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{d1} &= hf_1(x+h,y_1+k_{c1},y_2+k_{c2},\cdots,y_N+k_{cN}) \\
k_{d2} &= hf_2(x+h,y_1+k_{c1},y_2+k_{c2},\cdots,y_N+k_{cN}) \\
&\vdots \\
k_{dN} &= hf_N(x+h,y_1+k_{c1},y_2+k_{c2},\cdots,y_N+k_{cN}) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
最後に求めた\(k_a,k_b,k_c,k_d\)を使って\(x+h\)でのそれぞれの関数の値を導きます。、
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
x &= x+h \\
y_{1} &= y_{1}+(k_{a1}+2k_{b1}+2k_{c1}+k_{d1})/6 \\
y_{2} &= y_{2}+(k_{a2}+2k_{b2}+2k_{c2}+k_{d2})/6 \\
&\vdots \\
y_{N} &= y_{N}+(k_{aN}+2k_{bN}+2k_{cN}+k_{dN})/6 \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
となります。

4次ルンゲ・クッタ法のサブルーチン


fortran90でサブルーチンを書けばこんなコードになるかと思います。

subroutine rk4(N,x,h,y)
  !   sikinote!
  !N     --> Unknown number of first order differential equation
  !x     --> variable number of differential equations
  !y     --> function, differentiated by variable number(t)
  !
  !   /  dy1            
  !  |  ----- = f1(x,y1,y2,...)
  !  |    dx            
  !  |
  !  |   dy2            
  !  |  ----- = f2(x,y1,y2,...)
  !  |    dx            
  !  |
  !  |  .....
  !  |
  !  |   dyN            
  !  |  ----- = fN(x,y1,y2,...)
  !   \   dx            
  !
  implicit none
  integer::i,N
  double precision::h,rkfd,x
  double precision::tmp(1:N),y(1:N),K(1:4,1:N)
  external::rkfd

  tmp(1:N)=y(1:N)
  do i=1,N
     K(1,i)=h*rkfd(N,x,tmp,i)
  enddo

  x=x+0.5d0*h

  tmp(1:N)=y(1:N)+K(1,1:N)*0.5d0
  do i=1,N
     K(2,i)=h*rkfd(N,x,tmp,i)
  enddo

  tmp(1:N)=y(1:N)+K(2,1:N)*0.5d0
  do i=1,N
     K(3,i)=h*rkfd(N,x,tmp,i)
  enddo

  x=x+0.5d0*h

  tmp(1:N)=y(1:N)+K(3,1:N)
  do i=1,N
     K(4,i)=h*rkfd(N,x,tmp,i)
  enddo
 
  y(1:N)=y(1:N)+(K(1,1:N)+2.d0*K(2,1:N)+2.d0*K(3,1:N)+K(4,1:N))/6.d0
 
  return
end subroutine rk4

実際に例題を解きましょう。
2つの例題を解きます。


スポンサーリンク


1階微分方程式


1つ目は
\(
\displaystyle \frac{dy}{dx}=-y\sin{x}
\)

です。
一般解は
\(
y=Ae^{\cos{x}}
\)

であり、初期条件\(x=0,y=2\)として解けば解析解は
\(
y=2e^{\cos{x}-1}
\)

です。
これを計算するにはfortran90のプログラムを上の載せたサブルーチン”rk4″と共に、こう書けばいいです。

program main
  !sikinote
  implicit none
  integer::i,N,step,Nmax
  double precision::x,h
  double precision,allocatable::y(:)
 
  h=1.d-3
  Nmax=20000
  step=200
  N=1
  !############
  !# h    -> Step size of variable parameter
  !# Nmax -> Maximum number of calculation
  !# step -> Output interval of function
  !          (Number of output point -> Nmax/step)
  !# N    -> Number of first order differential equation
  !############
 
  allocate(y(1:N))
  !Input initial value to x,y(1:N)
  x=0.d0
  y(1)=2.d0
 
  !Output result to "y.txt"
  open(22,file='./y.txt',status='replace')
 
  !Output y(x=0)
  write(22,*)x,y(1)
  !Time propagation
  do i=1,Nmax
     call rk4(size(y),x,h,y)
     if(mod(i,step).eq.0)then
        write(22,*)x,y(1)
     endif
  enddo

  deallocate(y)
  close(22)
  stop
end program main
function rkfd(N,x,y,sign)
  !sikinote
  implicit none
  integer::N,sign
  double precision::rkfd,x,y(1:N)
 
  !Input right side function of first order differential equations.
  !This case correspond to below differential equations.
  !
  !   /  dy1            
  !  |  ----- = -y1*sin(x)
  !  |    dx            
  !   \
  !
  rkfd=0.d0
  if(sign.eq.1)then
     rkfd=-y(1)*sin(x)
  endif

  return
end function rkfd

解いて、gnuplotで
plot “y.txt”
数値計算解(赤点)解析解(緑線)と共に出力すると、こんなグラフが得られます。
4次ルンゲクッタ法による計算1

連立1階微分方程式


もう一つの例題は微分方程式
\(
\displaystyle \frac{d^2y}{dx^2}+2\gamma \frac{dy}{dx}+y=0
\)

を解きます。これは物理の世界ではバネの減衰振動を表す運動方程式で(詳しくは減衰振動へ。)、\(0\lt\gamma\lt1\)の場合、解は
\(
y(x)=Ae^{-\gamma x}\cos(x\sqrt{1-\gamma^2}-\alpha)
\)

となります。
\(\gamma=0.15\), 初期条件を今、
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\left.\frac{dy}{dx}\right|_{x=0} &= 1 \\
\left.y\right|_{x=0} &= -0.15 \\
\end{aligned}
\right.
\end{eqnarray}
\)
とすると、解析解は
\(
y(x)=e^{-\gamma x}\cos{(x\sqrt{1-0.15^2})}
\)

となります。
プログラムで計算する際は、まず連立1次微分方程式に焼きなおす必要があります。すなわち、
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{dy}{dx}&=v \\
\frac{dv}{dx}&=-2\gamma v -y \\
\end{aligned}
\right.
\end{eqnarray}
\)
として解けばいいわけです。

その時の変更するべき場所は少しで、mainとrkfdの一部を変えればおしまいです。
変えた場所をコメントで
!+—+
を入れたので確認してください。

program main
  !sikinote
  implicit none
  integer::i,N,step,Nmax
  double precision::x,h
  double precision,allocatable::y(:)
 
  h=1.d-3
  Nmax=20000
  step=200
  N=2           !+---+
  !############
  !# h    -> Step size of variable parameter
  !# Nmax -> Maximum number of calculation
  !# step -> Output interval of function
  !          (Number of output point -> Nmax/step)
  !# N    -> Number of first order differential equation
  !############
 
  allocate(y(1:N))
  !Input initial value to x,y(1:N)
  x=0.d0        !+---+
  y(1)=1.d0     !+---+
  y(2)=-0.15d0  !+---+
 
  !Output result to "y.txt"
  open(22,file='./y.txt',status='replace')
 
  !Output y(x=0)
  write(22,*)x,y(1),y(2)  !+---+
  !Time propagation
  do i=1,Nmax
     call rk4(size(y),x,h,y)
     if(mod(i,step).eq.0)then
        write(22,*)x,y(1),y(2)   !+---+
     endif
  enddo

  deallocate(y)
  close(22)
  stop
end program main
function rkfd(N,x,y,sign)
  !sikinote
  implicit none
  integer::N,sign
  double precision::rkfd,x,y(1:N)
 
  !Input right side function of first order differential equations.
  !This case correspond to below differential equations.
  !
  !   /  dy1            
  !  |  ----- = y2
  !  |    dx            
  !  |
  !  |   dy2            
  !  |  ----- = -0.3y2-y1
  !   \   dx      
  !
  rkfd=0.d0
  if(sign.eq.1)then
     rkfd=y(2)               !+---+
  elseif(sign.eq.2)then      !+---+
     rkfd=-0.3d0*y(2)-y(1)   !+---+
  endif

  return
end function rkfd

gnuplotで数値計算解(赤点)解析解(緑線)のグラフを書くとこうでなります。
4次ルンゲクッタ法による計算2

連立2階微分方程式


運動方程式に特化したサブルーチンです。
置いておきます。
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{d^2 x}{dt^2} &= -x-0.1v_x \\
\frac{d^2 y}{dt^2} &= -y-0.2v_y
\end{aligned}
\right.
\end{eqnarray}
\)

program main
  implicit none
  integer::i,N,step,Nmax
  double precision::t,h
  double precision,allocatable::r(:),v(:)
 
  h=1.d-3
  Nmax=20000
  step=200
  N=2
  !############
  !# h    -> Step size of variable parameter
  !# Nmax -> Maximum number of calculation
  !# step -> Output interval of function
  !          (Number of output point -> Nmax/step)
  !# N    -> Number of first order differential equation
  !############
 
  allocate(r(1:N),v(1:N))
  !Input initial value to t,r(1:N),v(1:N)
  t=0.d0
  r(1)=1.d0
  v(1)=-0.15d0
  r(2)=1.d0
  v(2)=0.5d0
 
  !Output result to "y.txt"
  open(22,file='./y.txt',status='replace')
 
  !Output y(x=0)
  write(22,*)t,r(1),r(2)
  !Time propagation
  do i=1,Nmax
     call rk4(size(r),t,h,r,v)
     if(mod(i,step).eq.0)then
        write(22,*)t,r(1),r(2)
     endif
  enddo

  deallocate(r)
  close(22)
  stop
end program main
function rkfd(N,t,r,v,sig,s)
  implicit none
  integer::N,sig,s
  double precision::rkfd
  double precision::t,r(1:N),v(1:N)
  !Input right side function of second order differential equations.
  !This case correspond to below differential equations.
  !
  !   /  d^2 x            
  !  |  ------- = -x-0.1v_x
  !  |    dt^2            
  !  |
  !  |   d^2 y            
  !  |  ------ = -y-0.2v_y
  !   \   dt^2      
  !
  rkfd=0.d0
  if(sig.eq.1)then
     rkfd=v(s)    
  elseif(sig.eq.2)then
     if(s.eq.1)then
        rkfd=-r(1)-0.1d0*v(1)
     elseif(s.eq.2)then
        rkfd=-r(2)-0.2d0*v(2)
     endif
  endif
 
  return
end function rkfd
subroutine rk4(N,t,h,r,v)
  !   sikinote!
  ! specialization for equation of motion.
  !N    --> Unknown number of second order differential equation
  !x    --> variable parameter of differential equations
  !r    --> function about position
  !v    --> function about velocity
  !
  !   /  d^2 y1            
  !  |  ------- = f1(x,y1,y2,...)
  !  |    dx^2            
  !  |
  !  |   d^2 y2            
  !  |  ------- = f2(x,y1,y2,...)
  !  |    dx^2            
  !  |
  !  |  .....
  !  |
  !  |   d^2 yN            
  !  |  ------- = fN(x,y1,y2,...)
  !   \   dx^2            
  !
  implicit none
  integer::i,N
  double precision::t,h
  double precision::rkfd,r(1:N),v(1:N),tmpr(1:N),tmpv(1:N),K(1:2,1:4,1:N)
  external::rkfd
 
  K=0.d0
  tmpr(1:N)=r(1:N)
  tmpv(1:N)=v(1:N)
  do i=1,N
     K(1,1,i)=h*rkfd(N,t,tmpr,tmpv,1,i)
     K(2,1,i)=h*rkfd(N,t,tmpr,tmpv,2,i)
  enddo
 
  t=t+0.5d0*h  

  tmpr(1:N)=r(1:N)+K(1,1,1:N)*0.5d0
  tmpv(1:N)=v(1:N)+K(2,1,1:N)*0.5d0
  do i=1,N
     K(1,2,i)=h*rkfd(N,t,tmpr,tmpv,1,i)
     K(2,2,i)=h*rkfd(N,t,tmpr,tmpv,2,i)
  enddo
 
  tmpr(1:N)=r(1:N)+K(1,2,1:N)*0.5d0
  tmpv(1:N)=v(1:N)+K(2,2,1:N)*0.5d0
  do i=1,N
     K(1,3,i)=h*rkfd(N,t,tmpr,tmpv,1,i)
     K(2,3,i)=h*rkfd(N,t,tmpr,tmpv,2,i)
  enddo
 
  t=t+0.5d0*h
 
  tmpr(1:N)=r(1:N)+K(1,3,1:N)
  tmpv(1:N)=v(1:N)+K(2,3,1:N)
  do i=1,N
     K(1,4,i)=h*rkfd(N,t,tmpr,tmpv,1,i)
     K(2,4,i)=h*rkfd(N,t,tmpr,tmpv,2,i)
  enddo
 
  r(1:N)=r(1:N)+(K(1,1,1:N)+2.d0*K(1,2,1:N)+2.d0*K(1,3,1:N)+K(1,4,1:N))/6.d0
  v(1:N)=v(1:N)+(K(2,1,1:N)+2.d0*K(2,2,1:N)+2.d0*K(2,3,1:N)+K(2,4,1:N))/6.d0
 
  return
end subroutine rk4
スポンサーリンク


コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です