📂수치해석

4차 룽게-쿠타 메소드

메소드 ¹

$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}}$ 지점의 데이터까지 사용한 것을 알 수 있다. 더 일반화된 설명은 같이보기를 참고하자.

직관적인 유도

위의 항들을 하나씩 풀어서 가장 간단한 오일러 메소드 $y_{n+1} = y_{n} + h f ( x_{n} , y_{n} ) $ 를 통해 생각해보자. (편의상 $\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) $$ 역시 백워드 오일러 메소드 $y_{n+1} = y_{n} + h f ( x_{n+1} , y_{n+1} )$ 에 의해 $\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] $$ 은 한 스텝 $h$ 를 진행시키기 위해 $x_{n+{{1} \over {2}} }$ 지점의 데이터까지 사용하면서 특히 $x_{n+{{1} \over {2}} }$ 지점에 가중치를 준 가중평균으로 볼 수 있다.

여기서 더 많은 계산을 추가할 순 있지만 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 플랏을 얻는다.

줄리아

다음은 위와 같은 내용의 코드를 줄리아로 구현한 것이다.

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)

같이보기

• 명시적 룽게-쿠타 메소드
• 암시적 룽게-쿠타 메소드