In [3]:
%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 [4]:
def lorenz_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 [5]:
assert np.allclose(lorenz_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 [6]:
def solve_lorenz(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*max_time))
soln=odeint(lorenz_derivs,ic,t,args=(sigma,rho,beta))
return soln,t
In [7]:
# solve_lorenz([5.0,0.0,1.0],max_time=0.1)
In [8]:
assert True # leave this to grade solve_lorenz
In [ ]:
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 [9]:
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 [10]:
# np.random.seed(1)
# icl=[np.random.randint(-15,15,10),np.random.randint(-15,15,10),np.random.randint(-15,15,10)]
# print(icl)
# solnsy=np.zeros((10,1))
# solnst=np.zeros((10,1))
# print('solntst:\n',solnst)
# l=list(icl[0])
# print(l[1])
# icx=list(icl[0])
# icy=list(icl[1])
# icz=list(icl[2])
# ic=[icx[1],icy[1],icz[1]]
# print(ic)
In [18]:
# def plot_lorenz(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)
# icl=[np.random.randint(-15,15,N),np.random.randint(-15,15,N),np.random.randint(-15,15,N)]
# print(icl)
# colors = plt.cm.hot(np.linspace(0,1,N))
# plt.figure(figsize=(11,8))
# icx=list(icl[0])
# icy=list(icl[1])
# icz=list(icl[2])
# for i in range(N):
# ic=[icx[i],icy[i],icz[i]]
# solnsy,solnst=solve_lorenz(ic,max_time,sigma,rho,beta)
# plt.plot(solnsy[0],solnsy[2],label='%s,%s '%(ic[0],ic[2]), color=colors[i])
# plt.legend(loc=4,title='Initial conditions(x,z)')
# plt.title('Trajectories')
# plt.ylim(-15.5,15.5)
# return solnsy
#the above plots lines but I dont know why I am going to try below with solving my lorenz first with all the values
#but this seems to waste power by storing all of them
#worked with Jessi Pilgram for the below attempt
def plot_lorenz(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)
icl=[np.random.randint(-15,15,N),np.random.randint(-15,15,N),np.random.randint(-15,15,N)]
colors = plt.cm.hot(np.linspace(0,1,N))
plt.figure(figsize=(11,8))
icx=list(icl[0])
icy=list(icl[1])
icz=list(icl[2])
solns=[]
for i in range(N):
ic=[icx[i],icy[i],icz[i]]
solnsy,solnst=solve_lorenz(ic,max_time,sigma,rho,beta)
solns.append(solnsy)
for i in solns:
solnx=[a[0] for a in i]
solnz=[b[2] for b in i]
plt.plot(solnx,solnz)
plt.legend(loc=4,title='Initial conditions(x,z)')
plt.title('Trajectories')
return solnsy
# # I was able to get colormap to work in
In [19]:
a=plot_lorenz()
In [20]:
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 [21]:
interact(plot_lorenz,max_time=(1,10,1),N=(1,50,1),sigma=(0.0,50.0,.1),rho=(0.0,50.0,.1),beta=fixed(8/3))
Out[21]:
Describe the different behaviors you observe as you vary the parameters $\sigma$, $\rho$ and $\beta$ of the system:
As N is increased the number of trajectories increases, which makes sense becasue the number of inital conditions are increasing. As max time increases the internal spirals increased. As sigma increases the spacing of the internal spirals increases but entire shape of the systme remains the same. As rho increases the internal rings increase while the entire systme flattens out.