Fourth-order Runge-Kutta method 📂Numerical Analysis

Fourth-order Runge-Kutta method

Method ¹

Given the initial value problem $\begin{cases} y ' = f(x,y) \\ y( x_{0} ) = Y_{0} \end{cases}$ for the continuous function $f$ defined in $D \subset \mathbb{R}^2$, let’s assume the interval $(a,b)$ is divided into nodes as in $a \le x_{0} < x_{1} < \cdots < x_{n} < \cdots x_{N} \le b$. Particularly, for sufficiently small $h > 0$, if $x_{j} = x_{0} + j h$, then for the initial value $y_{0} = Y_{0}$ y_{n+1} = y_{n} + h \sum_{j=1}^{p} \gamma_{j} V_{j} y_{n+1} = y_{n-1} + h \sum_{j=1}^{p} \gamma_{j} V_{j} $$

Explanation

The Runge-Kutta method, like the Adams method, comes in various forms and determines $\gamma_{j}$ and $V_{j}$ through complex algebraic operations. Among them, the 4th-order Runge-Kutta method, often abbreviated as RK4, is the most popularly used.

4th-order Runge-Kutta method: $$ \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*} $$

From the formula, you can see that to proceed one step $h$, data up to the point $\displaystyle {{h} \over {2}}$ is used. For a more generalized explanation, refer to See also.

Intuitive Derivation

Let’s break down the terms one by one and think using the simplest Euler method (Assuming $\displaystyle x_{n+ {{1} \over {2}} } : = x_{n} + {{h} \over {2}}$, $y_{n+ {{1} \over {2}} } \approx Y(x_{n + {{h} \over {2}} } )$ for convenience.) $$ V_{1} = f(x_{n} , y_{n} ) $$ is the differential coefficient $y_{n} ' $ at the point $x_{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*} $$ is the differential coefficient $\displaystyle y_{n + {{1} \over {2}} } ' $ at the point $x_{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) $$ is also by the backward Euler method, $\displaystyle V_{3} = f \left( x_{n + {{1} \over {2}} } , y_{n + {{1} \over {2}} } \right) $, the differential coefficient $\displaystyle y_{n + {{1} \over {2}} } ' $ at the point $x_{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*} $$ Furthermore, by the backward Euler method, it is the differential coefficient $\displaystyle y_{n + 1} ' $ at the point $x_{n+ 1 }$. $$ {{h} \over {6}} \left[ V_{1} + 2 V_{2} + 2 V_{3} + V_{4} \right] $$ is especially a weighted average that uses data up to the point $x_{n+{{1} \over {2}} }$ and gives a weight to the point $x_{n+{{1} \over {2}} }$ to proceed one step $h$.

Although more calculations can be added, from RK4 onwards, the order of convergence doesn’t increase as much as the amount of computation anymore. Even if not increased, if the problem is not stiff, RK4 is often sufficient, and its coding is also comparatively easy, which is why it is widely used.

Implementation

R

The following is the code to implement RK4 in R and draw the Lorenz attractor in 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')

Executing the code yields the following 3D plot.

$20190226\_173942.png$

Julia

Below is the implementation of the same content in Julia.

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)