Exercise 9.2: Numerically solving differential equations

---

import numpy as np

import bokeh.io
import bokeh.plotting

bokeh.io.output_notebook()

Loading BokehJS ...

For this exercise, we will use Euler’s method to simulate a classic set of differential equations that can exhibit chaotic solutions known as the Lorenz attractor. As a warm-up example, we will simulate bacterial growth.

Bacterial growth can be modeled by the differential equation

\begin{align} \frac{\mathrm{d}n}{\mathrm{d}t} = k n, \end{align}

where \(n\) is the number of bacteria and \(k\) is the growth rate. The idea here is that the number of bacteria will grow faster the more bacteria we have, because there are more to divide. Analytically, we know the solution to this differential equation is

\begin{align} n(t) = n_0 \mathrm{e}^{kt}, \end{align}

i.e., exponential growth. But suppose we did not know how to derive that. We could simulate the differential equation. We do this by discretizing time. Instead of a derivative, we have a change in \(n\) over a change in time \(t\).

\begin{align} \frac{\mathrm{d}n}{\mathrm{d}t} \approx \frac{\Delta n}{\Delta t} = k n. \end{align}

Let’s say we know \(n\) and time zero, \(n(0)\). Then \(n\) at time \(t = \Delta t\) is

\begin{align} n(\Delta t) \approx n(0) + \Delta n = n(0) + \Delta t\, k n(0). \end{align}

More generally, we can write

\begin{align} \frac{\mathrm{d}n}{\mathrm{d}t} = f(n), \end{align}

and

\begin{align} n(t+\Delta t) \approx n(t) + \Delta t\,f(n). \end{align}

So, we can instruct Python to take our current value of \(n\), and then add \(\Delta t\) times \(f(n)\) to get our new \(n\) at a time just a bit later on, at \(t + \Delta t\). We repeat this over and over again to move forward in time. Let’s code that up!

# Specify parameter
k = 1

# Specify my little time step
delta_t = 0.01

# Make an array of time points, evenly spaced up to 10
t = np.arange(0, 10, delta_t)

# Make an array to store the number of bacteria
n = np.empty_like(t)

# Set the initial number of bacteria
n[0] = 100

# Write a for loop to keep updating n as time goes on
for i in range(1, len(t)):
    n[i] = n[i - 1] + delta_t * k * n[i - 1]

Ok! We just computed the time points and the number of bacteria, so we can just plot the result!

p = bokeh.plotting.figure(
    height=300,
    width=450,
    x_axis_label="time (units of 1/k)",
    y_axis_label="number of bacteria",
)

p.line(x=t, y=n, line_width=2)

bokeh.io.show(p)

And there is the famous exponential growth!

This time stepping method is called Euler’s method, and what we’re doing is called numerical solution of a differential equation.

a) Now it’s time to solve the Lorenz attractor. In this case, there are three differential equations.

\begin{align} &\frac{\mathrm{d}x}{\mathrm{d}t} = \sigma(y - x),\\[1em] &\frac{\mathrm{d}y}{\mathrm{d}t} = x(\rho - z) - y,\\[1em] &\frac{\mathrm{d}z}{\mathrm{d}t} = xy - \beta z. \end{align}

Your task in this exercise is to numerically solve these two differential equations together and then plot the result. Use the following parameter values, which were the ones originally used by Lorenz.

sigma = 10
beta = 8/3
rho = 28

Start at time \(t = 0\) with \(x = y = z= 1\).

Even though there are now three differential equations, the procedure is the same, you update each by adding \(\Delta t\) times the respective derivative. You should integrate these equations until time \(t = 60\). I would recommend a step size of 0.001.

b) Plot your solution showing \(x\) versus time.

c) Solve again, this time with \(x(0) = 1.0001\), just 0.01% different from your previous starting point. Plot this solution on the same plot as your did in part (b) so you can compare the trajectories for just slightly different initial conditions.

d) Now make three plots. Plot \(x\) vs \(y\), \(x\) vs \(x\), and \(y\) vs \(z\). Look at the result. Cool, right?

e) [Bonus] Euler’s method is probably the simplest way to solve differential equations, and is by no means the best. SciPy has an ODE solver, scipy.integrate.odeint() that uses the more sophisticated and robust methods for solving systems of ODEs. Read the documentation about how scipy.integrate.odeint() works and use it to solve the Lorenz attractor system of ODEs.

This last part is tough; I’m not giving you directions, and you are kind of on your own to read the documentation and figure it out.