© 2018 Suzy Beeler and Griffin Chure. This work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license

This exercise was generated from a Jupyter notebook. You can download the notebook here.

In this tutorial, we will cover how to write a simple numerical integrator using the Forward Euler method to examine the dynamics of exponential growth.

Developing simple ways to solve ordinary differential equations has long been an area of intense research. While deriving the analytical solution may be simple in some cases, it is often useful to solve them numerically, especially if slamming out the analytical solution will give you carpal tunnel.

While there are many ways to numerically integrate these equations, in this tutorial we will examine the Forward Euler method. Say we have an ordinary differential equation such as

$$\frac{dN}{dt} = r N(t) \tag{1}$$as would be the case for exponential growth, where $r$ is some growth-rate constant and $t$ is time. Rather than solving this analytically (although it is trivial), we can solve it numerically by starting at some given value of $N$, evaluating Equation (1) for a given time step $\Delta t$, and updating the new value of $N$ at this new time $t+ \Delta t$. We can phrase this mathematically as

$$N(t+ \Delta t) = N(t) + rN(t) \Delta t .\tag{2}$$Say our initial value ($N$ at $t=0$) is $N=10$ and $r=1$. We can take a time step $\Delta t=0.1$ and find that the change in value of $N$ is

$$\Delta N = rN\Delta t = 1. \tag{3}$$We can then compute the new value of $N$ at time $t+\Delta t$ as

$$N(t+\Delta t) = N(t) + \Delta N = 10 + 1 = 11.\tag{4}$$We can then take another step forward in time and repeat the process for as long as we would like. As the total time we'd like to integrate over becomes large, it becomes obvious why using a computer is a far more attractive approach than scribbling it by hand.

A major point to be wary of is the instability of this method. The error in this scales with the square of our step size. We must choose a sufficiently small step in time such that at most only one computable event must occur. For example, if we are integrating exponential growth of bacterial cells, we don't want to take time steps larger than a cell division! This requirement is known as the Courant-Friedrichs-Lewy condition and is important for many different time-marching computer simulations.

As is often the case, the best way to learn is to do. Let's give our digits some exercise and numerically integrate this exponential growth differential equation.

In [ ]:

```
# Import the necessary modules
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# For pretty plots
import seaborn as sns
rc={'lines.linewidth': 2, 'axes.labelsize': 14, 'axes.titlesize': 14, \
'xtick.labelsize' : 14, 'ytick.labelsize' : 14}
sns.set(rc=rc)
```

In [ ]:

```
# parameters for our ODE
N_0 = 1
r = 0.03 # min^-1
# parameters for our integration
dt = 0.1 # min
total_time = 120 # min
```

`N_t`

array being set to `N_0`

as specified above.

In [ ]:

```
# determine the number of steps that will be taken
num_steps = int(total_time/dt)
# initilize an array of length num_steps into which to store values of N
N_t = np.zeros(num_steps)
N_t[0] = N_0
```

`N_t`

array and filling in the values of `N`

as we go. At each time point $N(t) = N(t-\Delta t) + rN(t-\Delta t) \Delta t,$ where $t - \Delta t$ refers to the previous entry in the `N_t`

array.

In [ ]:

```
# numerically integrate by looping through N_t
for t in range(1,num_steps):
# first calculate dN, using pevious N_t entry
dN = N_t[t-1] * r * dt
# update current N_t entry
N_t[t] = N_t[t-1] + dN
```

In [ ]:

```
# make array of time values
times = np.arange(num_steps)*dt
# plot
plt.plot(times,N_t)
plt.xlabel("time (mins)")
plt.ylabel("N")
```

Out[ ]:

Text(0, 0.5, 'N')

In [ ]:

```
# compute the known solution
soln = N_0 * np.exp(r*times)
# plot both our integration and the known solution
plt.plot(times,N_t)
plt.plot(times,soln)
plt.xlabel("time (mins)")
plt.ylabel("N")
plt.legend(["numerical integration", "known solution"])
```

Out[ ]:

<matplotlib.legend.Legend at 0x1a22aa0b50>

`dt`

. As `dt`

increases, our numerical integration deviates more and more from the known solution, with our integration systematically underestimating the true values. This is because with exponential growth, the rate of growth is always increasing and taking too large of a time step fails to capture this increase.