📂数値解析

四次のルンゲ＝クッタ法

メソッド ¹

$D \subset \mathbb{R}^2$ で定義された連続関数 $f$ に対する初期値問題 $\begin{cases} y ' = f(x,y) \\ y( x_{0} ) = Y_{0} \end{cases}$ が与えられている。区間 $(a,b)$ を $a \le x_{0} < x_{1} < \cdots < x_{n} < \cdots x_{N} \le b$ のようなノードポイントで分けたとする。特に充分に小さい $h > 0$ に対して $x_{j} = x_{0} + j h$ とすると、初期値 $y_{0} = Y_{0}$ に対して $$ y_{n+1} = y_{n} + h \sum_{j=1}^{p} \gamma_{j} V_{j} $$

説明

ルンゲ＝クッタ法は、アダムス法のように、多様な形式を持ち、複雑な代数的操作を通じて $\gamma_{j}$ と $V_{j}$ を決定する。その中でも特に人気があるのが4次のルンゲ＝クッタ法で、通常 RK4 と略される。

4次のルンゲ＝クッタ法: $$ \begin{align*} y_{n+1} =& y_{n} + {{h} \over {6}} \left[ V_{1} + 2 V_{2} + 2 V_{3} + V_{4} \right] \\ V_{1} =& f(x_{n} , y_{n}) \\ V_{2} =& f \left( x_{n} + {{h} \over {2}} , y_{n} + {{h} \over {2}} V_{1} \right) \\ V_{3} =& f \left( x_{n} + {{h} \over {2}} , y_{n} + {{h} \over {2}} V_{2} \right) \\ V_{4} =& f \left( x_{n} + h , y_{n} + h V_{3} \right) \end{align*} $$

この式を見ると、ステップ $h$ を進めるために $\displaystyle {{h} \over {2}}$ 地点のデータまで使ったことがわかる。もっと一般化した説明は関連項目を参照してほしい。

直感的な導入

上記の項を一つずつ解いて、最も単純なオイラー法を使って考えてみよう。 (便宜上 $\displaystyle x_{n+ {{1} \over {2}} } : = x_{n} + {{h} \over {2}}$、$y_{n+ {{1} \over {2}} } \approx Y(x_{n + {{h} \over {2}} } )$ とする。) $$ V_{1} = f(x_{n} , y_{n} ) $$ は、$x_{n}$ 地点での微分係数 $y_{n} ' $ になる。 $$ \begin{align*} V_{2} =& f \left( x_{n} + {{h} \over {2}} , y_{n} + {{h} \over {2}} f \left( x_{n} , y_{n} \right) \right) \\ =& f \left( x_{n + {{1} \over {2}} } , y_{n + {{1} \over {2}} } \right) \end{align*} $$ は、$x_{n+{{1} \over {2}} }$ 地点での微分係数 $\displaystyle y_{n + {{1} \over {2}} } ' $ になる。 $$ V_{3} = f \left( x_{n} + {{h} \over {2}} , y_{n} + {{h} \over {2}} f \left( x_{n + {{1} \over {2}} } , y_{n + {{1} \over {2}} } \right) \right) $$ また、バックワードオイラー法により $\displaystyle V_{3} = f \left( x_{n + {{1} \over {2}} } , y_{n + {{1} \over {2}} } \right) $ で、$x_{n+{{1} \over {2}} }$ 地点での微分係数 $\displaystyle y_{n + {{1} \over {2}} } ' $ になる。 $$ \begin{align*} V_{4} =& f \left( x_{n} + h , y_{n} + h V_{3} \right) \\ =& f \left( x_{n} + {{h} \over {2}} + {{h} \over {2}} , y_{n} + {{h} \over {2}} V_{3} + {{h} \over {2}} V_{3} \right) \\ =& f \left( x_{n + {{h} \over {2}}} + {{h} \over {2}} , y_{n + {{h} \over {2}}} + {{h} \over {2}} V_{3} \right) \end{align*} $$ さらに、バックワードオイラー法により、$x_{n+ 1 }$ 地点での微分係数 $\displaystyle y_{n + 1} ' $ になる。 $$ {{h} \over {6}} \left[ V_{1} + 2 V_{2} + 2 V_{3} + V_{4} \right] $$ 特に、$x_{n+{{1} \over {2}} }$ 地点のデータまで使いながら、特に $x_{n+{{1} \over {2}} }$ 地点に重みをつけた加重平均と見ることができるステップ $h$ を進めるためのものと考えられる。

ここで更に計算を加えることは可能だが、RK4からは収束次数が計算量だけで上がってくれない。別にもっと上げなくても、問題がスティフでなければ、RK4だけで十分であり、コーディングも比較的簡単で広く使われている。

実装

R

以下は、RでRK4を実装し、ローレンツアトラクタを3Dで描画するコードである。

library(rgl)
 
RK4<-function(ODE,v,h=10^(-2)){
  V1 = ODE(v)
  V2 = ODE(v + (h/2)*V1)
  V3 = ODE(v + (h/2)*V2)
  V4 = ODE(v + h*V3)
  v = v + (h/6)*(V1 + 2*V2 + 2*V3 + V4)
  return(v)
}
 
lorenz<-function(v,rho=28,sigma=10,beta=8/3){
  dvdt = v
  dvdt[1] = sigma*(v[2]-v[1])
  dvdt[2] = v[1]*(rho-v[3])-v[2]
  dvdt[3] = v[1]*v[2] - beta*v[3]
  return(dvdt)
}
 
traj<-numeric(0)
v=c(1,1,1)
for(i in 1:10000){
  traj<-rbind(traj,v)
  v=RK4(lorenz,v)
}
 
plot3d(traj,type='l')
plot(traj[,1],traj[,3],type='l')

このコードを実行すると、次のような3Dプロットを得られる。

ジュリア

../../categories/줄리아 以下は、上記と同じ内容のコードをジュリアで実装したものである。

using Plots
function RK4(ODE::Function,v::Array{Float64,1},h=10^(-2))
    V1 = ODE(v)
    V2 = ODE(v .+ (h/2)*V1)
    V3 = ODE(v .+ (h/2)*V2)
    V4 = ODE(v .+ h*V3)
    return @. v + (h/6)*(V1 + 2*V2 + 2*V3 + V4)
end
function lorenz(v::Array{Float64,1}; ρ=28.,σ=10.,β=8/3)
  dvdt = deepcopy(v)
  dvdt[1] = σ*(v[2]-v[1])
  dvdt[2] = v[1]*(ρ-v[3])-v[2]
  dvdt[3] = v[1]*v[2] - β*v[3]
  return dvdt
end
function  lorenz_attracter(v::Array{Float64,1}, endtime = 10000)
  x,y,z = [v[1]], [v[2]], [v[3]]
  for t in 1:endtime
    newx,newy,newz = RK4(lorenz,[x[end], y[end], z[end]])
    push!(x, newx)
    push!(y, newy)
    push!(z, newz)
  end
  return x,y,z
end
result = lorenz_attracter([1.,1.,1.])
plot(result)

関連項目

• 明示的ルンゲ＝クッタ法
• 暗示的ルンゲ＝クッタ法