The fourth order Runge-Kutta (RK4)
Enwu Liu
2023-05-03
The Runge-Kutta method attempts to improve the Euler’s method. Each Runge–Kutta method is derived from an appropriate Taylor method (In fact, Euler’s method can also be seen derived from a simple Taylor method). These methods can be constructed for any order. By far, the most often used method is the classical fourth-order Runge-Kutta method.
Suppose we want to solve the the following ODE
\[\left\{\begin{matrix} y'=f(t,y)\\ y(t_0)=y_0 \end{matrix}\right. \] For RK4 we will use the following recursive method to find the final solution
\[y_{n+1}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)h\] where: \[y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)h=\text{weighted average slope}\] \[\begin{array}{l} k_1 = f(t_n , y_n ), \\ k_2 = f\left( t_n + \frac{h}{2}, y_n + \frac{h}{2} k_1 \right) , \\ k_3 =f\left( t_n + \frac{h}{2}, y_n + \frac{h}{2} k_2 \right) , \\ k_4 = f( t_n + h, y_n + h k_3 ). \end{array}\]Other methods: Gill’s Algorithm RK4, Fifth order Butcher’s method, Runge–Kutta–Fehlberg Method (RKF45), et al
Example: R code to solve \[{{dy(t)} \over {dt}} + 2y(t) = 0\quad {\rm{or }} \quad {{dy(t)} \over {dt}} = - 2y(t)\] From: https://lpsa.swarthmore.edu/NumInt/NumIntFourth.html, I modified the Matlab code.
# Solve y'(t)=-2y(t) with y0=3, 4th order Runge Kutta
# Initial condition gives solution at t=0.
RK4_1<-function(h,t0,t_m,y0){
t = seq(t0,t_m,by=h)
nstep=(t_m-t0)/h+1
ystar<-c((1:nstep)*0)
ystar[1] = y0
for(i in 1:(length(t)-1)){
k1 = -2*ystar[i] # Approx for y gives approx for deriv
y1 = ystar[i]+k1*h/2 # Intermediate value (using k1)
k2 = -2*y1 # Approx deriv at intermediate value.
y2 = ystar[i]+k2*h/2 # Intermediate value (using k2)
k3 = -2*y2 # Another approx deriv at intermediate value.
y3 = ystar[i]+k3*h # Endpoint value (using k3)
k4 = -2*y3 # Approx deriv at endpoint value.
ystar[i+1] = ystar[i] + (k1+2*k2+2*k3+k4)*h/6 # Approx soln
}
return(list('ystar'=ystar,'t'=t))
}
y1<-RK4_1(0.2,0,2,3)
y2<-RK4_1(0.8,0,2,3)
t = seq(0,2,by=0.2) # t goes from 0 to 2 seconds.
yexact = 3*exp(-2*t)
plot(y1$t,y1$ystar,xlab='time',ylab='y')
lines(t,yexact,type='l')
lines(y2$t,y2$ystar,lty=3)
Another example solve \(\frac{dy}{dt}=y-t^2+1\) and \(y(0)=0.5\)
#f(t,y)=y-t^2+1
f_ty<-function(t,y){
f_ty<-y-t^2+1
return(f_ty)
}
RK4_2<-function(h,t0,t_m,y0){
t = seq(t0,t_m,by=h)
ystar<-c((1:length(t))*0)
ystar[1] = y0
for(i in 1:(length(t)-1)){
k1 = h*f_ty(t[i],ystar[i])
k2 = h*f_ty(t[i]+0.5*h,ystar[i]+k1/2)
k3 = h*f_ty(t[i]+0.5*h,ystar[i]+k2/2)
k4 = h*f_ty(t[i]+h,ystar[i]+k3)
ystar[i+1] = ystar[i] + (k1+2*k2+2*k3+k4)/6
}
return(ystar)
}
y<-RK4_2(0.5,0,2,0.5)
y## [1] 0.500000 1.425130 2.639603 4.006819 5.301605The exact solution of is \(y=t^2+2t+1-\frac{1}{2}e^t\)
t = seq(0,2,by=0.5)
y_exact<-t^2+2*t+1-0.5*exp(t)
y_exact## [1] 0.500000 1.425639 2.640859 4.009155 5.305472Runge-Kutta for a system of differential equations
\[ \left\{\begin{array}{l} \frac{dy}{dx}=f(x,y(x),z(x)), & y(x_0)=y_0 \\ \frac{dz}{dx}=g(x, y(x), z(x)),& z(x_0)=z_0\\ \end{array} \right. \]
\[\begin{array}{l} k_1 = h · f(x_n, y_n, z_n)\\ l_1 = h · g(x_n, y_n, z_n)\\ k_2 = h · f(x_n + h/2, y_n + k_1/2, z_n + l_1/2)\\ l_2 = h · g(x_n + h/2, y_n + k_1/2, z_n + l_1/2)\\ k_3 = h · f(x_n + h/2, y_n + k_2/2, z_n + l_2/2)\\ l_3 = h · g(x_n + h/2, y_n + k_2/2, z_n + l_2/2)\\ k_4 = h · f(x_n + h, y_n + k_3, z_n + l_3)\\ l_4 = h · g(x_n + h, y_n + k_3, z_n + l_3)\\ k = 1/6 · (k_1 + 2 · k_2 + 2 · k_3 + k_4)\\ l = 1/6 · (l_1 + 2 · l_2 + 2 · l_3 + l_4)\\ x_n+1 = x_n + h\\ y_n+1 = y_n + k\\ z_n+1 = z_n + l\\ \end{array}\]Example
\[\left\{\begin{array}{l} \frac{dy}{dx} = z \\ \frac{dz}{dx} = 6y - z \end{array}\right. \]
With \(y(0)=3\) and \(z(0)=1\), for \(0\le t \le 1\)
R code for solving above system using RK4:
F_xyz<-function(x,y,z){
fxyz<-z
return(fxyz)
}
G_xyz<-function(x,y,z){
gxyz<-6*y-z
return(gxyz)
}
RK4_3<-function(h,t0,t_m,y0,z0){
t = seq(t0,t_m,by=h)
ystar<-c((1:length(t))*0)
zstar<-c((1:length(t))*0)
ystar[1] = y0
zstar[1] = z0
for(i in 1:(length(t)-1)){
k1 = h*F_xyz(t[i],ystar[i],zstar[i])
l1 = h*G_xyz(t[i],ystar[i],zstar[i])
k2 = h*F_xyz(t[i]+0.5*h,ystar[i]+k1/2,zstar[i]+l1/2)
l2 = h*G_xyz(t[i]+0.5*h,ystar[i]+k1/2,zstar[i]+l1/2)
k3 = h*F_xyz(t[i]+0.5*h,ystar[i]+k2/2,zstar[i]+l2/2)
l3 = h*G_xyz(t[i]+0.5*h,ystar[i]+k2/2,zstar[i]+l2/2)
k4 = h*F_xyz(t[i]+h,ystar[i]+k3,zstar[i]+l3)
l4 = h*G_xyz(t[i]+h,ystar[i]+k3,zstar[i]+l3)
ystar[i+1] = ystar[i] + (k1+2*k2+2*k3+k4)/6
zstar[i+1] = zstar[i] + (l1+2*l2+2*l3+l4)/6
}
return(list('ystar'=ystar,'zstar'=zstar))
}
RK4_3(0.1,0,1,3,1)## $ystar
## [1] 3.000000 3.183638 3.532476 4.050814 4.752267 5.659661 6.805469
## [8] 8.232750 9.996623 12.166270 14.827579
##
## $zstar
## [1] 1.000000 2.663088 4.320751 6.068622 7.998407 10.203527 12.784314
## [8] 15.853107 19.539561 23.996434 29.406157We also can use R deSolve package to solve the above ODE system
ODE_sys<-function(t,state,parameters){
with(as.list(c(state)),
{
dy <-z
dz <- 6*y-z
list(c(dy,dz))
})
}
### define parameters
state <- c(y=3,z=1)
times<-seq(0,1,by=0.1)
require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=ODE_sys))
out## time y z
## 1 0.0 3.000000 1.000000
## 2 0.1 3.183624 2.663155
## 3 0.2 3.532461 4.320864
## 4 0.3 4.050808 6.068766
## 5 0.4 4.752278 7.998580
## 6 0.5 5.659696 10.203737
## 7 0.6 6.805536 12.784573
## 8 0.7 8.232863 15.853438
## 9 0.8 9.996792 19.539990
## 10 0.9 12.166514 23.996997
## 11 1.0 14.827920 29.406900We see the two methods got exactly the same results.
Another example:
\[\left\{\begin{array}{l} \frac{dy}{dt} = -y+4x \\ \frac{dx}{dt} = -y+6x \end{array}\right. \] with With \(y(0)=-0.5\) and \(x(0)=0.5\), for \(0\le t \le 1\) for \(0\le t \le 0.7\)
References
1.https://www.cfm.brown.edu/people/dobrush/am33/Matlab/ch3/RK4.html
2.https://mathworld.wolfram.com/Runge-KuttaMethod.html
3.https://math.okstate.edu/people/yqwang/teaching/math4513_fall11/Notes/rungekutta.pdf