Lab 11: Runge Kutta Methods of Higher Orders

In this lab, we implement a few higher order Runge Kutta methods for solving linear ODEs of the first order. Remember that an ODE of the first order involves only the first derivatives.

There are three popular methods for solving linear ODEs of the first order.

  1. Euler Method
  2. Modified Euler (or RK-1) Method
  3. RK-4 Method

In this lab, we implement the most popular fourth-order Runge Kutta Method (RK-4). RK-4 is popular because it has been found that higher order methods such as RK-5 are not as efficient.

In this method, we define the general form as: $$ \begin{eqnarray} k_1 & = & \Delta{t}f(t_n, x(t_n)) \\ k_2 & = & \Delta{t}f(t_n+\alpha\Delta{t}, x(t_n)+\beta{k}_1) \\ k_3 & = & \Delta{t}f(t_n+\gamma\Delta{t}, x(t_n)+\delta{k}_1+\epsilon{k}_2) \\ k_4 & = & \Delta{t}f(t_n+\zeta\Delta{t}, x(t_n) + \eta{k}_1+\theta{k}_2+\iota{k}_3) \end{eqnarray} $$ and finally, $$ x(t_{n+1}) = x(t_n) + ak_1 + bk_2 + ck_3 $$

Thankfully, the most popular RK-4 Method does not use so much of the Greek alphabet but really is quite simple! $$ \begin{eqnarray} k_1 & = & \Delta{t}f(t_n, x(t_n)) \\ k_2 & = & \Delta{t}f(t_n+\frac{\Delta{t}}{2}, x(t_n)+\frac{k_1}{2}) \\ k_3 & = & \Delta{t}f(t_n+\frac{\Delta{t}}{2}, x(t_n)+\frac{k_2}{2}) \\ k_4 & = & \Delta{t}f(t_n+\Delta{t}, x(t_n)+k_3),\,\,\,\,\, \mbox{and}\\ x(t_{n+1}) & = & x(t_n) + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{eqnarray} $$

Implement RK-4 in the cell below

In [1]:
import numpy as np
###################### ADD YOUR RK-4 CODE BELOW
def ode_rk4(func, x0, t0, tn, dt=0.2) :
    
    return xtn

Boundary Value Problems

Boundary value problems specify the value at two different locations (or times). These are considered the extremes of the range over which we are interested in the solution. A classic example is the differential equation $$ -\frac{d^2x}{dt^2} = f(t) $$ with $x = A$ at $t=t_0$ and $x = B$ at $t=t_1$. The problem is to find values of $x$ between $t_0$ and $t_1$.

We saw that the solution is obtained by defining $N$ points between the two boundaries (including the boundaries). Then it is solved in our oh-so-loved way using gauss_eliminate() and back_substitute() methods!

The solution is done in three steps:

  1. Define the $A$ matrix as $$ \left(\begin{array}{ccccccc} 2 & 1 & 0 & 0 & \ldots & 0 & 0 \\ -1& 2 & -1& 0 & \ldots & 0 & 0 \\ 0 & -1& 2 & -1& \ldots & 0 & 0 \\ \vdots & & \vdots & & \ldots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \ldots & -1& 2 \end{array}\right) $$
  2. Define the array $Y$ as $$ \left(\begin{array}{c} h^2f(t_1) + A \\ h^2f(t_2) \\ h^2f(t_3) \\ \vdots \\ h^2f(t_{n-1}) \\ h^2f(t_n) + B \end{array}\right) $$
  3. Solve the equation $AX = Y$

Implement the boundary value problem in the cell(s) below.

In [2]:
def ode_bvalue(func, x1, x2, t1, t2, h=0.2) :
    
    return xvals