The Lorenz Equations

Consider the Lorenz equations:

formula not implemented

Import relevant Python packages

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
1.5s

Define a function that evaluates the right-hand side of the Lorenz equations

The function takes in two inputs: time and 'state'. 'State' is a vector of [x,y,z] values. It returns the rate of change of the state, which is given by the right-hand side of an equation of the form inline_formula not implemented. We leave 'time' as one of the input variables accepted by the function to keep things general; in this particular system, inline_formula not implemented has no time-dependence, and that input is not used by the function.

def lorenz_rhs(t,state,s,r,b):
    # Break up state variable into 3:
    x,y,z = state
    dxdt = s*(y-x)
    dydt = r*x - y - x*z
    dzdt = x*y - b*z
    # Assemble
    return [dxdt,dydt,dzdt]
0.0s

First, the usual values

Set up Lorenz's usual parameters

We choose an arbitrary initial condition; it doesn't matter much here, because basically any point in phase space will eventually end up on the 'strange attractor' discovered by Lorenz.

We choose a time-span of 50 time units, and plot it 10,000 times to get nice and smooth curves.

# Initial condition
y0 = [0.9, 0.0, 0.0]
# Time span for the solution (start time, end time)
t_span = (0, 50)
# Plot the solution at the following times:
tvals = np.linspace(0,t_span[1],10000)
# Solve the ODE
solution1 = solve_ivp(lorenz_rhs, # system of equations
                      t_span, # time span of integration
                      y0, # initial conditions
                      args=(10,28,8/3), # parameters
                      t_eval=tvals, # times at which to store solution
                      method='RK45' # method of integration
                     )
0.3s

Plot Lorenz's solution

On the left-hand side, we plot the full 3-dimensional phase space. On the right-hand side, we plot one of the variables against time.

fig1 = plt.figure(dpi=100,figsize=(12,5))
# 3-d plot
ax1_1 = fig1.add_subplot(1,2,1, projection='3d')
ax1_1.plot(solution1.y[0], solution1.y[1], solution1.y[2], lw=0.5, color='black')
ax1_1.set_title("Solution in [x,y,z] space ")
# time-trace
ax1_2 = fig1.add_subplot(1, 2, 2)
ax1_2.set_aspect(1/2)
ax1_2.plot(solution1.t, solution1.y[2], lw=0.5, color='blue')
ax1_2.set_title("Variable 'z' against time")
plt.show()
0.7s

Now for some new values: try r = 100

This time, create a new `solve_ivp` object, using different integration times and parameter values as needed.

Solve the ODE again

# Solve the ODE
t_end2 = 70
solution2 = solve_ivp(lorenz_rhs, # system of equations
                      (0, t_end2), # time span of integration
                      [0.0, 1.0, 0.0], # initial conditions
                      args=(10,100,8/3), # parameters
                      t_eval=np.linspace(0,t_end2,10000), # times at which to store solution
                     )
0.6s

Then, plot the solution

fig2 = plt.figure(dpi=100,figsize=(12,5))
# 3-d plot
ax2_1 = fig2.add_subplot(1,2,1, projection='3d')
ax2_1.plot(solution2.y[0], solution2.y[1], solution2.y[2], lw=0.5, color='black')
ax2_1.set_title("Solution in [x,y,z] space ")
# time-trace
ax2_2 = fig2.add_subplot(1, 2, 2)
ax2_2.set_aspect(0.1)
ax2_2.plot(solution2.t, solution2.y[2], lw=0.5, color='blue')
ax2_2.set_title("Variable 'z' against time")
plt.show()
0.8s

Take a closer look at the last few time units:

In Python, you can solve over a large range of times and choose to plot only a subset of those times. Let us use this fact to plot the last 1,000 time units only:

plt.plot(solution2.t[-1000:],solution2.y[2][-1000:],lw=0.5)
0.7s

What do you notice?

Next, try: r = 101, r = 150, r = 170, and r = 215

For these, focus on the time-trace plots of any of the three variables x,y, and/or z.

Runtimes (1)