```
In [1]:
```import numpy as np

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