Open Educational Resources

Solution Methods for IVPs: Runge-Kutta Methods

Runge-Kutta Methods

The Runge-Kutta methods developed by the German mathematicians C. Runge and M.W. Kutta are essentially a generalization of all the previous methods. Consider the following IVP:

    \[\frac{\mathrm{d}x}{\mathrm{d}t}=F(x,t)\]


Assuming that the value of the dependent variable x (say x_i) is known at an initial value t_i, then, the Runge-Kutta methods employ the following general equation to calculate the value of x at t=t_{i+1}, namely x_{i+1} with h=t_{i+1}-t_i:

    \[x_{i+1}=x_i+h\phi\]


where \phi is a function of the form:

    \[\phi=\alpha_1 k_1+\alpha_2k_2+\alpha_3k_3+\cdots+\alpha_nk_n\]


\alpha_j are constants while k_j=F(x,t) evaluated at points within the interval [x_i,x_{i+1}] and have the form:

    \[\begin{split}k_1&=F(x_i,t_i)\\k_2&=F(x_i+q_{11}k_1h,t_i+p_1h)\\k_3&=F(x_i+q_{21}k_1h+q_{22}k_2h,t_i+p_2h)\\&\vdots\\k_n&=F(x_i+q_{(n-1),1}k_1h+q_{(n-1),2}k_2h+\cdots+q_{(n-1),(n-1)}k_{n-1}h,t_i+p_{n-1}h)\end{split}\]

The constants \alpha_j, and the forms of k_j are obtained by equating the value of x_{i+1} obtained using the Runge-Kutta equation to a particular form of the Taylor series. The k's are recurrence relationships, meaning k_1 appears in k_2, which appears in k_3, and so forth. This makes the method efficient for computer calculations. The error in a particular form depends on how many terms are used. The general forms of these Runge-Kutta methods could be implicit or explicit. For example, the explicit Euler method is a Runge-Kutta method with only one term with \alpha_1=1 and k_1=F(x_i,t_i). Heun’s method on the other hand is a Runge-Kutta method with the following non-zero terms:

    \[\begin{split}\alpha_1&=\alpha_2=\frac{1}{2}\\k_1&=F(x_i,t_i)\\k_2&=F(x_i+hk_1,t_i+h)\end{split}\]


Similarly, the midpoint method is a Runge-Kutta method with the following non-zero terms:

    \[\begin{split}\alpha_2&=1\\k_1&=F(x_i,t_i)\\k_2&=F(x_i+\frac{h}{2}k_1,t_i+\frac{h}{2})\end{split}\]


The most popular Runge-Kutta method is referred to as the “classical Runge-Kutta method”, or the “fourth-order Runge-Kutta method”. It has the following form:

    \[x_{i+1}=x_i+\frac{h}{6}\left(k_1+2k_2+2k_3+k_4\right)\]


with

    \[\begin{split}k_1&=F(x_i,t_i)\\k_2&=F\left(x_i+\frac{h}{2}k_1,t_i+\frac{h}{2}\right)\\k_3&=F\left(x_i+\frac{h}{2}k_2,t_i+\frac{h}{2}\right)\\k_4&=F(x_i+hk_3,t_i+h)\end{split}\]


The error in the fourth-order Runge-Kutta method is similar to that of the Simpson’s 1/3 rule, with a local error term E_{\mbox{local}}=\mathcal{O}(h^5) which would translate into a global error term of E=\mathcal{O}(h^4). The method is termed fourth-order because the error term is directly proportional to the step size raised to the power of 4.

View Mathematica Code
RK4Method[fp_, x0_, h_, t0_, tmax_] := (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[k1 = fp[xtable[[i - 1]], t0 + (i - 2) h];
   k2 = fp[xtable[[i - 1]] + h/2*k1, t0 + (i - 1.5) h];
   k3 = fp[xtable[[i - 1]] + h/2*k2, t0 + (i - 1.5) h];
   k4 = fp[xtable[[i - 1]] + h*k3, t0 + (i - 1) h];
   xtable[[i]] = xtable[[i - 1]] + h (k1 + 2 k2 + 2 k3 + k4)/6, {i, 2,
     n}];
  Data2 = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data2)

View Python Code
def RK4Method(fp, x0, h, t0, tmax):
  n = int((tmax - t0)/h + 1)
  xtable = [0 for i in range(n)]
  xtable[0] = x0
  for i in range(1,n):
    k1 = fp(xtable[i - 1], t0 + (i - 1)*h)
    k2 = fp(xtable[i - 1] + h/2*k1, t0 + (i - 0.5)*h)
    k3 = fp(xtable[i - 1] + h/2*k2, t0 + (i - 0.5)*h)
    k4 = fp(xtable[i - 1] + h*k3, t0 + i*h)
    xtable[i] = xtable[i - 1] + h*(k1 + 2*k2 + 2*k3 + k4)/6 
  Data = [[t0 + i*h, xtable[i]] for i in range(n)]
  return Data

Example

Solve Example 4 above using the classical Runge-Kutta method but with h=0.4.

Solution

Recall the differential equation of this problem:

    \[\frac{\mathrm{d}x}{\mathrm{d}t}=F(x,t)=tx^2+2x\]


Setting x_0=-5, t_0=0, h=0.4, the following are the values of k_1, k_2, k_3, and k_4 required to calculate x_1:

    \[\begin{split}k_1&=F(x_0,t_0)=F(-5,0)\\&=0+2(-5)=-10\\k_2&=F\left(x_0+\frac{h}{2}k_1,t_0+\frac{h}{2}\right)=F(-5+0.2(-10),0+0.2)\\&=0.2(-5+0.2(-10))^2+2\times (-5+0.2(-10))=-4.2\\k_3&=F\left(x_0+\frac{h}{2}k_2,t_0+\frac{h}{2}\right)=F(-5+0.2(-4.2),0+0.2)\\&=0.2(-5+0.2(-4.2))^2+2\times (-5+0.2(-4.2))=-4.8589\\k_4&=F(x_0+hk_3,t_0+h)=F(-5+0.4(-4.8589),0+0.4)\\&=0.4(-5+0.4(-4.8589))^2+2\times(-5+0.4(-4.8589))=5.3981\end{split}\]


Therefore:

    \[x_1=x_0+\frac{h}{6}(k_1+2k_2+2k_3+k_4)=-5+\frac{0.4}{6}(-10-2\times 4.2-2\times 4.8589+5.3981)=-6.51465\]

Proceeding iteratively gives the values of x_i up to t=5.2. The Microsoft Excel file RK4.xlsx provides the required calculations. The following Microsoft Excel table shows the required calculations:

RK42

The following Mathematica code implements the classical Runge-Kutta method for this problem with h=0.3. The output curve is shown underneath.

View Mathematica Code
RK4Method[fp_, x0_, h_, t0_, tmax_] := (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[k1 = fp[xtable[[i - 1]], t0 + (i - 2) h];
   k2 = fp[xtable[[i - 1]] + h/2*k1, t0 + (i - 1.5) h];
   k3 = fp[xtable[[i - 1]] + h/2*k2, t0 + (i - 1.5) h];
   k4 = fp[xtable[[i - 1]] + h*k3, t0 + (i - 1) h];
   xtable[[i]] = xtable[[i - 1]] + h (k1 + 2 k2 + 2 k3 + k4)/6, {i, 2,
     n}];
  Data2 = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data2)
Clear[x, xtable]
a = DSolve[{x'[t] == t*x[t]^2 + 2*x[t], x[0] == -5}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := t*x^2 + 2*x;
Data2 = RK4Method[fp, -5.0, 0.3, 0, 5];
Plot[x, {t, 0, 5}, Epilog -> {PointSize[Large], Point[Data2]},  AxesOrigin -> {0, 0}, AxesLabel -> {"t", "x"}, PlotRange -> All]
Title = {"t_i", "x_i"};
Data2 = Prepend[Data2, Title];
Data2 // MatrixForm

View Python Code
# UPGRADE: need Sympy 1.2 or later, upgrade by running: "!pip install sympy --upgrade" in a code cell
# !pip install sympy --upgrade
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
sp.init_printing(use_latex=True)

def RK4Method(fp, x0, h, t0, tmax):
  n = int((tmax - t0)/h + 1)
  xtable = [0 for i in range(n)]
  xtable[0] = x0
  for i in range(1,n):
    k1 = fp(xtable[i - 1], t0 + (i - 1)*h)
    k2 = fp(xtable[i - 1] + h/2*k1, t0 + (i - 0.5)*h)
    k3 = fp(xtable[i - 1] + h/2*k2, t0 + (i - 0.5)*h)
    k4 = fp(xtable[i - 1] + h*k3, t0 + i*h)
    xtable[i] = xtable[i - 1] + h*(k1 + 2*k2 + 2*k3 + k4)/6 
  Data = [[t0 + i*h, xtable[i]] for i in range(n)]
  return Data

x = sp.Function('x')
t = sp.symbols('t')
sol = sp.dsolve(x(t).diff(t) - t*x(t)**2 - 2*x(t), ics={x(0): -5})
display(sol)
def fp(x, t): return t*x**2 + 2*x
Data = RK4Method(fp, -5.0, 0.3, 0, 5)
x_val = np.arange(0,5,0.01)
plt.plot(x_val, [sol.subs(t, i).rhs for i in x_val])
plt.scatter([point[0] for point in Data],[point[1] for point in Data],c='k')
plt.xlabel("t"); plt.ylabel("x")
plt.grid(); plt.show()
print(["t_i", "x_i"],"\n",np.vstack(Data))

RK41

The following tool provides a comparison between the explicit Euler method and the classical Runge-Kutta method. The classical Runge-Kutta method gives excellent predictions even when h=0.5, at which point the explicit Euler method fails to predict anything close to the exact equation. The classical Runge-Kutta method also provides better estimates than the midpoint method and Heun’s method above.




Lecture Video

 

 


 

 

 


 

Leave a Reply

Your email address will not be published. Required fields are marked *