#### Time grid
We start by discretising time, that is dividing our total time interval $[0, t_f]$ into $N$ sub-intervals of equal size.
Specifically, we consider a time grid $t_0 = 0, t_1, t_2, \cdots, t_n, \cdots, t_N = t_f$, where $t_n = n\Delta t$. $\Delta t$ is known as the **time step**. Note that there are $N+1$ grid points in time corresponding to $N$ time steps of size $\Delta t$. Many authors use $h$ for the time step, but we will stick to $\Delta t$.
In our Python code we will use
- `dt` for the time step $\Delta t$,
- `Nsteps` for the number of time steps $N$.
#### Euler time stepping
You will agree that
$
\frac{dy}{dt}(t) \simeq \frac{y(t+\Delta t) - y(t)}{\Delta t}
$
We use this to approximate the derivative $\displaystyle{\frac{dy}{dt}}$ on the left-hand-side of our [[Numerical Solutions of ODEs|ODE]]
$
\frac{y(t+\Delta t) - y(t)}{\Delta t} \simeq f(y(t),t)
$
and solving for $y(t+\Delta t)$
$
y(t+\Delta t) \simeq y(t) + \Delta t f(y(t),t)
$
*Euler time stepping* is obtained by replacing the true solution $y(t_n)$ by the numerical approximation $y_n$, and dropping the approximately equals sign
$
y_{n+1} = y_n + \Delta t f(y_n, t_n) \tag{2}
$
This is a rule for going from $y_n$ to $y_{n+1}$, a process known as taking one time step. Euler time stepping is in fact the simplest possible rule consistent with the ODE.
We assume that we are told the initial condition. Hence, we know $y_0$. We can simply repeatedly apply the Euler time stepping to find our approximate solution $y_n$ at later times.
# Basic NumPy Implementation
### Example 1: Exponential Growth ODE
Below is a complete Python code for solving the linear, one-variable ODE
$
\dot y = -\mu y
$
by Euler time stepping. $\mu$ is a parameter. The numerical solution and true solution are plotted. You can read and easily understand how it works given the previous discussion. Two comments:
- We use variable `y0` for the initial condition. After the solution array `y` is created, we initialise this array with `y[0]=y0`. This helps keep the problem specification, including setting the initial condition, separate from the construction of the numerical solution.
- The variable `ydot` might appear to be unnecessary, but it is useful to keep the ODE specific parts of the code separate from the numerical-method specific parts of the code (the Euler time step). Later, we will compute `ydot` using functions.
```run-python
# import libraries
import numpy as np
import matplotlib.pyplot as plt
# Solve the ODE: dot y = -mu y, by Euler time stepping,
# and plot the numerical and true solutions.
# Problem setup.
# Set model and numerical parameters, and the initial condition.
# These are the lines the user might want to vary.
tf = 5
Nsteps = 50
mu = 0.75
y0 = 5
# Generate the time grid and solution array
t, dt = np.linspace(0, tf, Nsteps+1, retstep=True) # retstep allows linspace method to return the time step
y = np.zeros(Nsteps+1)
# Set the first point in the solution array to the initial condition
y[0] = y0
# Euler time-stepping loop:
for n in range(Nsteps):
# ydot = RHS evaluated at t_n
ydot = -mu * y[n]
y[n+1] = y[n] + dt * ydot
# plot the numerical solution
plt.plot(t, y, 'o', label='numerical')
# plot the true solution
y_true = y0*np.exp(-mu*t)
plt.plot(t, y_true, label='true')
# labels etc
plt.xlabel("t", fontsize=14)
plt.ylabel("y", fontsize=14)
plt.legend()
plt.title("Euler solution to $\dot y = -\mu y
quot;, fontsize=16)
plt.show()
```
### Example 2: Damped Harmonic Oscillator
Consider the damped harmonic oscillator ODEs
$
\begin{align}
\dot y_1 & = y_2 \\
\dot y_2 & = - y_1 - \mu y_2,
\end{align}
$
where $\mu \ge 0$ is a parameter. For $\mu = 0$ these equations become the un-damped harmonic oscillator.
```run-python
# import libraries import numpy as np import matplotlib.pyplot as plt
# Set problem parameters, including IC and mu
Nsteps = 1000
tf = 6*np.pi
y0 = np.array([5, -1])
mu = 1
# Compute number of dependent variables from y0
Neqs = len(y0)
# Allocate t and y_sol arrays
t, dt = np.linspace(0, tf, Nsteps+1, retstep=True)
y_sol = np.zeros(shape=(Nsteps+1, Neqs))
# Initialise the solution array
y_sol[0,:] = y0[:]
def damped(y, t):
# returns RHS of damped harmonic oscillator ODE
y1, y2 = y
mu = 0.25
y1dot = y2
y2dot = -y1 - mu * y2
return y1dot, y2dot
# Euler time-stepping loop (for the identity function f(t,y) = y)
for n in range(Nsteps):
ydot = np.array(damped(y_sol[n,:],dt*n))
y_sol[n+1,:] = y_sol[n,:] + dt * ydot
# plot numerical solution
y1 = y_sol[:,0]
y2 = y_sol[:,1]
plt.plot(t,y1,'-', label="y1")
plt.plot(t,y2,'-', label="y2")
# labels etc
plt.xlabel("t", fontsize=14)
plt.ylabel("y1, y2", fontsize=14)
plt.legend(fontsize=14)
plt.title("Damped harmonic oscillator", fontsize=16)
plt.show()
# plot a phase portrait with markers at beginning and end.
plt.figure(figsize=(5,5))
plt.plot(y1, y2)
plt.plot(y1[0], y2[0],'d', markersize = 10)
plt.plot(y1[-1], y2[-1],'o', markersize = 10)
plt.xlim(-6, 6)
plt.ylim(-6, 6)
plt.xlabel("y1", fontsize=14)
plt.ylabel("y2", fontsize=14)
plt.title("Phase portrait", fontsize=16)
plt.show()
```
# SciPy Implementation using odeint
### Example 1: Damped Harmonic Oscillator
```run-python
# import libraries
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
def damped(y, t):
# returns RHS of damped harmonic oscillator ODE
y1, y2 = y
mu = 0.25
y1dot = y2
y2dot = -y1 - mu * y2
return y1dot, y2dot
# set initial conditions
y0 = np.array([5, -1])
# set up time grid for solution
Nsteps = 1000
tf = 6 * np.pi
t= np.linspace(0, tf, Nsteps+1)
# We do not create the y_sol array. odeint will do that for us.
# call solver
y_sol = odeint(damped, y0, t)
y1 = y_sol[:,0]
y2 = y_sol[:,1]
# plot the numerical solution
plt.plot(t,y1,'-', label="y1")
plt.plot(t,y2,'-', label="y2")
# labels etc
plt.xlabel("t", fontsize=14)
plt.ylabel("y1, y2", fontsize=14)
plt.legend(fontsize=14)
plt.title("Damped harmonic oscillator", fontsize=16)
plt.show()
# plot a phase portrait with markers at beginning and end.
plt.figure(figsize=(5,5))
plt.plot(y1, y2)
plt.plot(y1[0], y2[0],'d', markersize = 10)
plt.plot(y1[-1], y2[-1],'o', markersize = 10)
plt.xlim(-6, 6)
plt.ylim(-6, 6)
plt.xlabel("y1", fontsize=14)
plt.ylabel("y2", fontsize=14)
plt.title("Phase portrait", fontsize=16)
plt.show()
```
### Example 2: [[SEIR Model]]