

In [1]:

    
import numpy as np

Steppers

Various steppers which are used to numerically solve differential equations are derived and defined in this notebook

The function below numerically solves an equation of the form $$\frac{dy}{dt}=f(y,t)$$

by approximating using the Euler method

$$ \begin{align*} y_{n+1}&=y_{n}+h \dfrac{dy}{dt}\\ y_{n+1}&=y_{n}+h f(y_n,t_n) \end{align*} $$

Where $y$ is the state vector of the system

The derivation is as follows



In [2]:

    
#Solves an equation of the form dy/dt=f(y,t) using the Euler method
#f is the function which returns the derivative of the state vector of the system
#y0 is the initial state of the system
#t is an array of times at which f will be evaluated

def Euler(f, y0, t, args=()):
    n    = len(t)                 #Find number of steps
    y0   = np.asarray(y0)         #Convert input state to numpy array
    y    = np.zeros((n, len(y0))) #Initialise array to hold state at each timestep
    y[0] = y0                     #Set initial state of system
        
    for i in range(n-1):
        h      = t[i+1]-t[i]
        y[i+1] = y[i] + h * f(y[i], t[i], *args) #*args passes all members of args to function
    return y

The function below numerically solves an equation of the form $$\frac{dy}{dt}=f(y,t)$$

by approximating using the Backward Euler method.

$$ \begin{align*} y_{n+1}&=y_{n}+h f(y_{n+1},t_{n+1}) \end{align*} $$

This approximation yields an algebraic equation for $y_{n+1}$ in which $y_{n+1}$ appears on both sides. This is because the Backward Euler method is an implicit method. The implicit equation must be solved by fixed point iteration (as used in this case) or some version of the Newton-Rhapson method.

The the function $f(y,t)$ must be evaluated multiple times each step with implicit methods. However implicit methods are also more stable than explicit methods.



In [3]:

    
#Solves an equation of the form dy/dt=f(y,t) using the backward Euler method
#f is the function which returns the derivative of the state vector of the system
#y0 is the initial state of the system
#t is an array of times at which f will be evaluated

def BackwardEuler(f, y0, t, args=()):
    n    = len(t)                 #Find number of steps
    y0   = np.asarray(y0)         #Convert input state to numpy array
    y    = np.zeros((n, len(y0))) #Initialise array to hold state at each timestep
    y[0] = y0                     #Set initial state of system
        
    for i in range(n-1):
        h      = t[i+1]-t[i]
        #Fixed point iteration to approximate y[i+1]
        y[i+1] = y[i]
        for j in range(10): 
            y[i+1] = y[i] + h * f(y[i+1], t[i+1], *args)  #*args passes all members of args to function
    return y

The function below numerically solves an equation of the form $$\frac{dy}{dt}=f(y,t)$$

by approximating using the improved Euler method (Heun's method) $$ \begin{align*} y_{n+1}&=y_{n}+\frac{h}{2}[f(y_n, t_n) + f(y_{n+1},t_{n+1})] \end{align*} $$

This method leads to a cancelation of errors as the forward and backward approximations of the derivative are averaged



In [4]:

    
#Solves an equation of the form dy/dt=f(y,t) using the improved Euler method, also know as Heun's Method
#f is the function which returns the derivative of the state vector of the system
#y0 is the initial state of the system
#t is an array of times at which f will be evaluated

def ImprovedEuler(f, y0, t, args=()):
    n    = len(t)                 #Find number of steps
    y0   = np.asarray(y0)         #Convert input state to numpy array
    y    = np.zeros((n, len(y0))) #Initialise array to hold state at each timestep
    y[0] = y0                     #Set initial state of system
        
    for i in range(n-1):
        h      = t[i+1]-t[i]
        k1     = f(y[i],     t[i],    *args) #Forward evaluation
        k2     = f(y[i]+h*k1,t[i]+h,  *args) #Backward evaluation
        y[i+1] = y[i] + h*(k1+k2)/2          #Central evaluation 
    return y

The function below numerically solves an equation of the form $$\frac{dy}{dt}=f(y,t)$$

by approximating using the Runge–Kutta method $$ y_{n+1} = y_{n}+\frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) $$ Where $$ \begin{align*} k_1 &= f(y_n, t_n)\\ k_2 &= f(y_n + \frac{h}{2}k_1, t_n + \frac{h}{2})\\ k_3 &= f(y_n + \frac{h}{2}k_2, t_n + \frac{h}{2})\\ k_4 &= f(y_n + h k_3, t_n + h)\\ \end{align*} $$

This method is a fourth order method so the local truncation error is of order $h^5$



In [5]:

    
#Solves an equation of the form dy/dt=f(y,t) using a Runge-Kutta method
#f is the function which returns the derivative of the state vector of the system
#y0 is the initial state of the system
#t is an array of times at which f will be evaluated

def RK4(f, y0, t, args=()):
    n    = len(t)                 #Find number of steps
    y0   = np.asarray(y0)         #Convert input state to numpy array
    y    = np.zeros((n, len(y0))) #Initialise array to hold state at each timestep
    y[0] = y0                     #Set initial state of system
        
    for i in range(n-1):
        h  = t[i+1] - t[i]
        k1 = f(y[i],  t[i],             *args)
        k2 = f(y[i] + h*k1/2, t[i]+h/2, *args)
        k3 = f(y[i] + h*k2/2, t[i]+h/2, *args)
        k4 = f(y[i] + h*k3,   t[i+1],   *args)
        y[i+1] = y[i] + h*(k1 + (2*k2) + (2*k3) + k4)/6
    return y