Lab 10: Ordinary Differential Equations

Ordinary Differential Equations (ODE) have functions and their derivatives as terms in the equations. For example, $$ \frac{d^2x}{dt^2} + 2t\frac{dx}{dt} + 4 = 0 $$ is an ordinary differential equation of the second order as it contains a second derivative. Note that we use $x$ and $t$ as the variables because many ODEs deal with functions and values that evolve over time. The above equation describes how the variable $x$ evolves over time. ODEs occur everywhere in scientific applications and many analytical methods have been developed for solving them. Equally importantly, it has been recognised that there are large families of ODEs which cannot be solved by analytical methods. The only option is numerical methods.

There are three broad sets of problems involving ODEs.

  1. Initial Value problems
  2. Boundary Value problems
  3. Eigenvalue problems
In this lab, we address the first, i.e. Initial value problems. In these equations, the value of the function is known at a single point, called the initial value. For example, the value of $x$ may be known at time $t_0$ as $x_0$ and it is needed to find the value after several time steps at time $t = t_n$. It is conventional to assume that the time step is a constant given by $\Delta{t}$ so that $t_n = t_{n-1} + \Delta{t} = t_0 + (n-1)\Delta{t}$.

Euler Method

Euler method is the simplest numerical method for solving ODEs. It applies to first order ODEs of the form $$ \frac{dx}{dt} = f(x, t). $$ It is fairly easy to solve the above with an approximation (derived from Taylor series expansions) $$ \begin{eqnarray} x(t_{n+1}) & = & x(t_n) + \left.\Delta{t} \frac{dx}{dt}\right|_{x = x(t_n)}\\ & = & x(t_n) + \Delta{t}f(x(t_n), t_n) \end{eqnarray} $$ As the value $x_0$ is known at $t_0$, the iterative equation above provides a solution.

Implement Euler method using the skeleton in the cell below.

In [ ]:
import numpy as np
from matplotlib import pyplot as plt
In [ ]:
def ode_euler(func, x0, t0, tn, dt=0.2) :
    
    return xt

Test your code on the following problem: Find $x(1)$ for $$ \frac{dx}{dt} = x $$ with $x(t_0) = 1$ at $t_0 = 0$. Try it initially with $dt = 0.2$. Later, reduce dt = 0.02 and dt = 1e-3. In the last case, $x(1) \approx e$ ($e = 2.71828182\ldots$).

In [ ]:
### Try your code here on the test problem given

Modified Euler or Runge-Kutta Order 1 Method

Euler method may be modified using a very simple trick to make it significantly better. Instead of using the derivative at $x = x(t_n)$, substitute it with the average of the slopes at $x = x(t_n)$ and $x = x(t_{n+1})$. $$ x(t_{n+1}) = x(t_n) + \Delta{t}\frac{f(x(t_n), t_n) + f(x(t_{n+1}), t_n + \Delta{t})}{2} $$ The problem, of course, is the $x(t_{n+1})$ term on the right-side which is unknown and precisely what we are aiming to compute.

Again, the trick is to estimate or predict the value using Euler method and then updating it with the Modified Euler equation in a two-step process. It is also called a predictor-corrector process. $$ \begin{eqnarray} x^p(t_{n+1}) & = & x(t_n) + \Delta{t}f(x(t_n), t_n) \,\,\,\,\,\mbox{(Predictor Step)} \\ k_1 & = & \Delta{t}f(x(t_n), t_n) \\ k_2 & = & \Delta{t}f(x^p(t_{n+1}), t_n + \Delta{t}) \\ & = & \Delta{t}f(x(t_n) + k_1, t_n + \Delta{t}) \\ x(t_{n+1}) & = & x(t_n) + \frac{1}{2}\left(k_1 +k_2\right) \,\,\,\,\,\,\,\mbox{(Corrector Step)} \end{eqnarray} $$ When we look at the higher order Runge-Kutta methods, we understand why this form of the modified Euler method is a first-order Runge-Kutta method.

Implement the Runge-Kutta order 1 method in the cell below.

In [ ]:
def ode_rk1(func, x0, t0, tn, dt=0.2) :
    
    return xtn

Do the same problem as above with the Euler method to verify your code. You will find that the method gives you a much better approximation to $e$.

In [ ]:
def func(x, t) :
    return x
In [ ]:
xinit = 1.0
tinit = 0.0
delta = 0.2
tfinal = 1.0
x1 = ode_euler(func, xinit, tinit, tfinal, dt=delta)
print(x1)
In [ ]:
xrk = ode_rk1(func, xinit, tinit, tfinal, dt=delta)
print(xrk)