In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from IPython.html.widgets import interact, fixed
The Lorenz system is one of the earliest studied examples of a system of differential equations that exhibits chaotic behavior, such as bifurcations, attractors, and sensitive dependence on initial conditions. The differential equations read:
$$ \frac{dx}{dt} = \sigma(y-x) $$$$ \frac{dy}{dt} = x(\rho-z) - y $$$$ \frac{dz}{dt} = xy - \beta z $$The solution vector is $[x(t),y(t),z(t)]$ and $\sigma$, $\rho$, and $\beta$ are parameters that govern the behavior of the solutions.
Write a function lorenz_derivs that works with scipy.integrate.odeint and computes the derivatives for this system.
In [2]:
def lorentz_derivs(yvec, t, sigma, rho, beta):
"""Compute the the derivatives for the Lorentz system at yvec(t)."""
x=yvec[0]
y=yvec[1]
z=yvec[2]
dx=sigma*(y-x)
dy=x*(rho-z)-y
dz=x*y-beta*z
return np.array([dx,dy,dz])
In [3]:
assert np.allclose(lorentz_derivs((1,1,1),0, 1.0, 1.0, 2.0),[0.0,-1.0,-1.0])
Write a function solve_lorenz that solves the Lorenz system above for a particular initial condition $[x(0),y(0),z(0)]$. Your function should return a tuple of the solution array and time array.
In [4]:
def solve_lorentz(ic, max_time=4.0, sigma=10.0, rho=28.0, beta=8.0/3.0):
"""Solve the Lorenz system for a single initial condition.
Parameters
----------
ic : array, list, tuple
Initial conditions [x,y,z].
max_time: float
The max time to use. Integrate with 250 points per time unit.
sigma, rho, beta: float
Parameters of the differential equation.
Returns
-------
soln : np.ndarray
The array of the solution. Each row will be the solution vector at that time.
t : np.ndarray
The array of time points used.
"""
t=np.linspace(0,max_time,int(250.0*max_time))
soln=odeint(lorentz_derivs,ic,t,args=(sigma,rho,beta))
return soln,t
In [5]:
assert True # leave this to grade solve_lorenz
Write a function plot_lorentz that:
N different initial conditions. To generate your initial conditions, draw uniform random samples for x, y and z in the range $[-15,15]$. Call np.random.seed(1) a single time at the top of your function to use the same seed each time.hot colormap from Matplotlib.The following cell shows how to generate colors that can be used for the lines:
In [6]:
N = 5
colors = plt.cm.hot(np.linspace(0,1,N))
for i in range(N):
# To use these colors with plt.plot, pass them as the color argument
print(colors[i])
In [7]:
def plot_lorentz(N=10, max_time=4.0, sigma=10.0, rho=28.0, beta=8.0/3.0):
"""Plot [x(t),z(t)] for the Lorenz system.
Parameters
----------
N : int
Number of initial conditions and trajectories to plot.
max_time: float
Maximum time to use.
sigma, rho, beta: float
Parameters of the differential equation.
"""
np.random.seed(1)
ic=np.random.rand(N,3)*30-15
plt.figure(figsize=(9,6))
# This takes the solutions of solve_lorentz in the x and z position of the
# array and uses the initial conditions of their respective positions.
for i in ic:
plt.plot(solve_lorentz(i,max_time,sigma,rho,beta)[0][:,0],solve_lorentz(i,max_time,sigma,rho,beta)[0][:,2]);
# I could not find a way to make the color mapping work
plt.xlabel('x(t)'),plt.ylabel('z(t)');
plt.title('Lorentz Parametric System')
In [8]:
plot_lorentz();
In [9]:
assert True # leave this to grade the plot_lorenz function
Use interact to explore your plot_lorenz function with:
max_time an integer slider over the interval $[1,10]$.N an integer slider over the interval $[1,50]$.sigma a float slider over the interval $[0.0,50.0]$.rho a float slider over the interval $[0.0,50.0]$.beta fixed at a value of $8/3$.
In [10]:
interact(plot_lorentz, max_time=[1,10], N=[1,50], sigma=[0.0,50.0], rho=[0.0,50.0], beta=fixed(8/3));
Describe the different behaviors you observe as you vary the parameters $\sigma$, $\rho$ and $\beta$ of the system:
$\bullet$ As $\sigma$ is varied, the shapes of the motion either become larger or smaller and move closer together as $\sigma$ increases.
$\bullet$ As $\rho$ is increased, the trajectory of each path circles more and more around a point on the graph.
$\bullet$ $\beta$ is fixed.
In [ ]: