In [18]:
%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 [19]:
def lorentz_derivs(yvec, t, sigma, rho, beta):
"""Compute the the derivatives for the Lorentz system at yvec(t)."""
# YOUR CODE HERE
#raise NotImplementedError()
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 [20]:
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 [21]:
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.
"""
# YOUR CODE HERE
#raise NotImplementedError()
t = np.linspace(0,max_time, max_time*250)
sol = odeint(lorentz_derivs, ic, t, args=(sigma, rho, beta))
return np.array(sol),np.array(t)
In [22]:
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 [23]:
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 [32]:
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.
"""
# YOUR CODE HERE
#raise NotImplementedError()
np.random.seed(1)
sols=[]
for i in range(N):
data = (np.random.random(3)-0.5)*30
sols.append(solve_lorentz(data, max_time, sigma, rho, beta))
for s in sols:
x = [p[0] for p in s[0]]
z = [p[2] for p in s[0]]
color = plt.cm.hot((x[0]+z[0])/60+0.5)
plt.plot(x,z, color=color)
plt.xlabel('$x(t)$')
plt.ylabel('$z(t)$')
plt.title('Lorentz System')
plt.box(False)
plt.tick_params(axis='x', top='off')
plt.tick_params(axis='y', right='off')
plt.xlim(-28,28)
plt.ylim(-20,65)
In [33]:
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]:
# YOUR CODE HERE
#raise NotImplementedError()
interact(plot_lorentz, max_time=(1,10,1), N = (1,50,1), 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:
When $\rho$ increases the trajectory expands in both x and z directions. It esentially tightens or loosens the orbital trajectory. $\sigma$ slightly expands the graph. It also, once under about 4.5, changes the side on which the main orbit occurs when it hits zero the trajectory is just a vericle line. $\beta$ is a fixed variable for this assigment. $N$ changes the number of trajectories plotted and max_time plots more or less points.