<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <table>
 <tr align=left><td><img align=left src="./images/CC-BY.png">
 <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Marc Spiegelman,  Based on ipython notebook by Kyle Mandli from his course [Introduction to numerical methods](https
Step1:  Exploring the Lorenz Equations
The Lorenz Equations are a 3-D dynamical system that is a simplified model of Rayleigh-Benard thermal convection.  They are derived and described in detail in Edward Lorenz' 1963 paper Deterministic Nonperiodic Flow  in the Journal of Atmospheric Science.  In their classical form they can be written
$$ 
    \dot{X} = \sigma( Y - X)\
    \dot{Y} = rX - Y - XZ \
    \dot{Z} = XY -b Z
$$
where $\sigma$ is the "Prandtl Number",  $r = \mathrm{Ra}/\mathrm{Ra}_c$ is a scaled "Raleigh number" and $b$ is a parameter that is related to the the aspect ratio of a convecting cell in the original derivation.
This ipython notebook, will provide some simple python routines for numerical integration and visualization of the Lorenz Equations.  The primary code is modified from the Wikipedia page
Some useful python routines
Step2:  Subcritical behavior $r<1$
Here we will begin exploring the behavior of the Lorenz equations for fixed values of $\sigma$ and $b$ and just changing the Rayleigh number $r$.
We will begin with subcritical behavior $r=0.5$  which rapidly damps to a condition of no motion
Step3:  Damped Oscillation $r=10$
Now we increase the Rayleigh number to  $r=10$ which admits two steady solutions depending on initial condition.
Step4:  Chaos and the strange attractor $r=28$
Now we increase the Rayleigh number to  $r=28$ and the solution becomes highly time-dependent and a-periodic.
Step5:  Limit Cycle at large Rayleigh number
Now we increase the Rayleigh number to  $r=350$ and the solution goes to a periodic limit cycle
Step6:  Stability of Fixed Points
It is straightforward to show that the Lorenz system has a fixed point at $X=Y=Z=0$ for all values of parameters $r,\sigma,b$.  Moreover, the Jacobian for the origin is
$$ 
    J = \left[
        \begin{matrix}
            -\sigma & \sigma & 0 \
                 r & -1 & 0      \ 
                 0 & 0 & -b 
        \end{matrix}
        \right]
$$
which has three negative eigenvalues for $r<1$ (Stable sink),  and 2 negative and 1 positive eigenvalue for $r > 1$ (3-D Saddle point).
At $r=1$, the Jacobian is singular and the origin undergoes a pitchfork bifurcation where two new fixed points appear at coordinates 
$$
    C^{\pm} = \left[
            \begin{matrix} 
            \pm\sqrt{b(r-1)} \
            \pm\sqrt{b(r-1)} \
            (r-1)
            \end{matrix}
            \right] 
$$
The stability of $C^+$ and $C^-$ depend on the eigenvalues of the Jacobian at these points, e.g 
$$ 
    J = \left[
        \begin{matrix}
            -\sigma & \sigma & 0 \
                 1 & -1 & -x^+      \ 
                 x^+ & x^+ & -b 
        \end{matrix}
        \right]
$$
where $x^+=\sqrt{b(r-1)}$.  These eigenvalues can be found as the roots of the cubic polynomial in $\lambda$ given by $|J -\lambda I|=0$ (which I will ask you to find in a homework problem).  But here we will just calculate and plot them numerically
Step7:  Liapunov Exponents
for $r>r_H$,  the Lorenz system exhibits a chaotic or "strange" attractor with nonperiodic orbits that exhibit extreme sensitivity to initial conditions.  One measure of this sensitivity is the "Liapunov Exponent" $\lambda$ defined roughly by
$$
        ||\delta(t)|| = ||\delta_0||e^{\lambda t}
$$ 
where $||\delta(t)||$ is the distance between to trajectories that started with initial conditions an infinitesimal distance $\delta_0<<1$ apart at $t=0$.
Here we will just estimate numerically the Liaponuv exponent for our chaotic system.
Step8:  Estimating the Liapunov exponent
Now we will calculate and plot the distance between the two trajectories on a log plot
Step9:  Calculating the Lorenz Map
The Lorenz Map is a discrete dynamical system that attemps to predict the $n+1$ maximum $Z_{n+1}(t)$ as a function of the previous peak $Z_{n}$.  Viewing the dynamical system in either the $Y,Z$ plane or as the time series $Z(t)$, we see that there are roughly two types of behavior,  a set of growing oscillations with increasing amplitude in $Z$, separated by jumps to a smaller value of $Z$ where it appears qualitatively, that during the jumps, the larger the value of $Z_n$, the smaller the value $Z_{n+1}$ after the jump.
Step10:  The Lorenz Map
The Lorenz map is simply a plot of the peak amplitude of the $n+1$ maximum $Z_{n+1}(t)$ as a function of the previous peak $Z_{n}$. Here we extract the maxima with a bit of discrete hackery where we
approximate the derivative of Z'
approximate extrema using intervals of the discrete solution with zero crossings of Z'
approximate local maxima using extrema with Z greater than mean(Z)
given the maxima $Z_n$,  plot $Z_{n+1}$ vs $Z_{n}$.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
from numpy.linalg import eigvals
Explanation: <table>
 <tr align=left><td><img align=left src="./images/CC-BY.png">
 <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Marc Spiegelman,  Based on ipython notebook by Kyle Mandli from his course [Introduction to numerical methods](https://github.com/mandli/intro-numerical-methods)</td>
</table>
End of explanation
def Lorenz(state,t,sigma,r,b):
  '''
  Returns the RHS of the Lorenz equations
  '''
  # unpack the state vector
  x = state[0]
  y = state[1]
  z = state[2]
  # compute state derivatives
  xd = sigma * (y-x)
  yd = (r-z)*x - y
  zd = x*y - b*z
  # return the state derivatives
  return [xd, yd, zd]
def SolveLorenz(state0,t,sigma=10.,r=28.,b=8./3.0):
    '''
    use ODEINT to integrate the lorenz equations from initial condition state0 at t=0 for
    the range of times given in the numpy array t
    '''
    Lorenz_p = lambda state,t: Lorenz(state,t,sigma,r,b)
    state = odeint(Lorenz_p, state0, t)
    return state
def PlotLorenzXvT(state,t,sigma,r,b):
    '''
    make time series plots of solutions of the Lorenz equations X(t),Y(t),Z(t)
    '''
    plt.figure()
    ax = plt.subplot(111)
    X = state[:,0]
    Y = state[:,1]
    Z = state[:,2]
    ax.plot(t,X,'r',label='X')
    ax.hold(True)
    ax.plot(t,Y,'g',label='Y')
    ax.plot(t,Z,'b',label='Z')
    ax.set_xlabel('time t')
    plt.title('Lorenz Equations: $\sigma=${}, $r=${}, $b=${}'.format(sigma,r,b))
    # Shrink current axis's height by 10% on the bottom
    box = ax.get_position()
    ax.set_position([box.x0, box.y0 + box.height * 0.1,
                 box.width, box.height * 0.9])
    # Put a legend below current axis
    ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.05),ncol=3)
    plt.show()
    
  
    
def PlotLorenz3D(state,sigma,r,b):
    '''
    Show 3-D Phase portrait using mplot3D
    '''
    # do some fancy 3D plotting
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot(state[:,0],state[:,1],state[:,2])
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    plt.title('Lorenz Equations: $\sigma=${}, $r=${}, $b=${}'.format(sigma,r,b))
    plt.show()
    
Explanation: Exploring the Lorenz Equations
The Lorenz Equations are a 3-D dynamical system that is a simplified model of Rayleigh-Benard thermal convection.  They are derived and described in detail in Edward Lorenz' 1963 paper Deterministic Nonperiodic Flow  in the Journal of Atmospheric Science.  In their classical form they can be written
$$ 
    \dot{X} = \sigma( Y - X)\
    \dot{Y} = rX - Y - XZ \
    \dot{Z} = XY -b Z
$$
where $\sigma$ is the "Prandtl Number",  $r = \mathrm{Ra}/\mathrm{Ra}_c$ is a scaled "Raleigh number" and $b$ is a parameter that is related to the the aspect ratio of a convecting cell in the original derivation.
This ipython notebook, will provide some simple python routines for numerical integration and visualization of the Lorenz Equations.  The primary code is modified from the Wikipedia page
Some useful python routines
End of explanation
# Set the parameters
sigma= 10.
b = 8./3
# set the initial condition
X0 = [2.0, 3.0, 4.0]
# set the time for integration
t = np.arange(0.0, 30.0, 0.01)
# set the Rayleigh number
r = 0.5
# solve the Equations
state = SolveLorenz(X0,t,sigma,r,b)
# and Visualize as a time series
PlotLorenzXvT(state,t,sigma,r,b)
# and as a 3-D phase portrait
PlotLorenz3D(state,sigma,r,b)
Explanation: Subcritical behavior $r<1$
Here we will begin exploring the behavior of the Lorenz equations for fixed values of $\sigma$ and $b$ and just changing the Rayleigh number $r$.
We will begin with subcritical behavior $r=0.5$  which rapidly damps to a condition of no motion
End of explanation
# set the Rayleigh number
r = 10.0
X0 = [2.,3.,4.]
state = SolveLorenz(X0,t,sigma,r,b)
PlotLorenzXvT(state,t,sigma,r,b)
PlotLorenz3D(state,sigma,r,b)
# now change the initial condition so X=-2
X0 = [-2.0, -3.0, 4.0]
state = SolveLorenz(X0,t,sigma,r,b)
PlotLorenzXvT(state,t,sigma,r,b)
PlotLorenz3D(state,sigma,r,b)
Explanation: Damped Oscillation $r=10$
Now we increase the Rayleigh number to  $r=10$ which admits two steady solutions depending on initial condition.
End of explanation
# set the Rayleigh number
r = 28.0
X0 = [2.,3.,4.]
state = SolveLorenz(X0,t,sigma,r,b)
PlotLorenzXvT(state,t,sigma,r,b)
PlotLorenz3D(state,sigma,r,b)
Explanation: Chaos and the strange attractor $r=28$
Now we increase the Rayleigh number to  $r=28$ and the solution becomes highly time-dependent and a-periodic.
End of explanation
# set the Rayleigh number
r = 350
X0 = [2.,3.,4.]
t = np.arange(0,8.,.0001)
state = SolveLorenz(X0,t,sigma,r,b)
PlotLorenzXvT(state,t,sigma,r,b)
PlotLorenz3D(state,sigma,r,b)
Explanation: Limit Cycle at large Rayleigh number
Now we increase the Rayleigh number to  $r=350$ and the solution goes to a periodic limit cycle
End of explanation
sigma = 10
b = 8./3.
r_H = sigma*(sigma+b+3)/(sigma-b -1.) # critical value of r at Hopf bifurcation
r_max = 28.
ra = np.linspace(1,28.,20)
xstar = lambda r: np.sqrt(b*(r-1))
J = lambda r: np.array([[-sigma,sigma,0],[1,-1,-xstar(r)],[xstar(r),xstar(r),-b]])
# plot out the eigenvalues
import  matplotlib.cm as cm
cmap = cm.get_cmap('coolwarm')
fig = plt.figure()
for r in ra:
    L = eigvals(J(r))
    plt.plot(np.real(L),np.imag(L),'o',color=cmap((r-min(ra))/(max(ra)-min(ra))))
    plt.hold(True)
# plot out eigenvalues at the Hopf Bifurcation
L = eigvals(J(r_H))
plt.plot(np.real(L),np.imag(L),'sy')
plt.xlabel('Re$(\lambda)$')
plt.ylabel('Im$(\lambda)$')
plt.title('Eigenvalues of $C^+$ for $r\in[1,{}]$, $r_H={}$'.format(max(ra),r_H))
plt.grid()
plt.show()
Explanation: Stability of Fixed Points
It is straightforward to show that the Lorenz system has a fixed point at $X=Y=Z=0$ for all values of parameters $r,\sigma,b$.  Moreover, the Jacobian for the origin is
$$ 
    J = \left[
        \begin{matrix}
            -\sigma & \sigma & 0 \
                 r & -1 & 0      \ 
                 0 & 0 & -b 
        \end{matrix}
        \right]
$$
which has three negative eigenvalues for $r<1$ (Stable sink),  and 2 negative and 1 positive eigenvalue for $r > 1$ (3-D Saddle point).
At $r=1$, the Jacobian is singular and the origin undergoes a pitchfork bifurcation where two new fixed points appear at coordinates 
$$
    C^{\pm} = \left[
            \begin{matrix} 
            \pm\sqrt{b(r-1)} \
            \pm\sqrt{b(r-1)} \
            (r-1)
            \end{matrix}
            \right] 
$$
The stability of $C^+$ and $C^-$ depend on the eigenvalues of the Jacobian at these points, e.g 
$$ 
    J = \left[
        \begin{matrix}
            -\sigma & \sigma & 0 \
                 1 & -1 & -x^+      \ 
                 x^+ & x^+ & -b 
        \end{matrix}
        \right]
$$
where $x^+=\sqrt{b(r-1)}$.  These eigenvalues can be found as the roots of the cubic polynomial in $\lambda$ given by $|J -\lambda I|=0$ (which I will ask you to find in a homework problem).  But here we will just calculate and plot them numerically
End of explanation
# start by running the Lorenz system long enough to get on the attractor
r = 28.0
X0 = [1.,0.,0.]
t = np.arange(0,20,.01)
state = SolveLorenz(X0,t,sigma,r,b)
# extract the final state and perturb it by a small amount epsilon
X0 = state[-1]
epsilon=1.e-6
X1 = X0 + epsilon*np.random.rand(3)
delta_0 = np.sqrt(np.sum((X1-X0)**2))
# Now run both initial conditions
t=np.arange(0.,50.,.0001)
state0 = SolveLorenz(X0,t,sigma,r,b)
state1 = SolveLorenz(X1,t,sigma,r,b)
# Compare the two trajectories as time-series X
plt.figure()
ax = plt.subplot(111)
ax.plot(t,state0[:,0],'r',t,state1[:,0],'b')
plt.xlabel('t')
plt.ylabel('X(t)')
plt.show()
# and in the phase space
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(state0[:,0],state0[:,1],state0[:,2],'r-')
plt.hold(True)
ax.plot(state1[:,0],state1[:,1],state1[:,2],'b-')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.title('Lorenz Equations: $\sigma=${}, $r=${}, $b=${}'.format(sigma,r,b))
plt.show()
Explanation: Liapunov Exponents
for $r>r_H$,  the Lorenz system exhibits a chaotic or "strange" attractor with nonperiodic orbits that exhibit extreme sensitivity to initial conditions.  One measure of this sensitivity is the "Liapunov Exponent" $\lambda$ defined roughly by
$$
        ||\delta(t)|| = ||\delta_0||e^{\lambda t}
$$ 
where $||\delta(t)||$ is the distance between to trajectories that started with initial conditions an infinitesimal distance $\delta_0<<1$ apart at $t=0$.
Here we will just estimate numerically the Liaponuv exponent for our chaotic system.
End of explanation
# calculate the distance between the two trajectories
delta = state1-state0
delta = np.sqrt(np.sum(delta**2,1))
# and plot them
plt.figure()
plt.semilogy(t,delta)
plt.xlabel('t')
plt.ylabel('$||\delta(t)||$')
plt.grid()
# now fit the first part with a straight line to determine the slope
# we'll pick the line between tmin and tmax to avoid initial transients and later saturation
tmin = 1.
tmax = 12.
imin = int(np.argwhere(t<tmin)[-1])
imax = int(np.argwhere(t>tmax)[0])
tfit = t[imin:imax]
p= np.polyfit(tfit,np.log(delta[imin:imax]),1)
plt.hold(True)
plt.semilogy(tfit,np.exp(p[1]+p[0]*tfit),'r')
plt.title('Liapunov Exponent Estimate $\lambda={}$'.format(p[0]))
plt.show()
Explanation: Estimating the Liapunov exponent
Now we will calculate and plot the distance between the two trajectories on a log plot
End of explanation
Y = state0[:,1]
Z = state0[:,2]
plt.figure()
plt.subplot(2,1,1)
ihalf = int(len(Y)/2.)
print ihalf,len(Y)
plt.plot(Y[:ihalf],Z[:ihalf])
plt.xlabel('Y')
plt.ylabel('Z')
plt.title('Lorenz system, $Y,Z$ plane: $r={}$, $\sigma={}$,$b={}$'.format(r,sigma,b))
plt.grid()
plt.subplot(2,1,2)
plt.plot(t,Z)
plt.xlabel('$t$')
plt.ylabel('$Z(t)$')
plt.title('$Z$ time series')
plt.show()
Explanation: Calculating the Lorenz Map
The Lorenz Map is a discrete dynamical system that attemps to predict the $n+1$ maximum $Z_{n+1}(t)$ as a function of the previous peak $Z_{n}$.  Viewing the dynamical system in either the $Y,Z$ plane or as the time series $Z(t)$, we see that there are roughly two types of behavior,  a set of growing oscillations with increasing amplitude in $Z$, separated by jumps to a smaller value of $Z$ where it appears qualitatively, that during the jumps, the larger the value of $Z_n$, the smaller the value $Z_{n+1}$ after the jump.
End of explanation
# first let's estimate the centered derivatve of Z to isolate the extrema 
dZ = np.zeros(Z.shape)
dZ[1:-2] = Z[2:-1] - Z[0:-3]
dZ.shape
plt.figure()
plt.plot(t,dZ,t,np.zeros(t.shape),'k:')
plt.ylabel('$dZ$')
plt.xlabel('t')
# now let's find all all intervals that contain zero crossings
icross = np.nonzero(dZ[:-2]*dZ[1:-1] <= 0)
Zextreme = Z[icross]
# and pick out all Extremes greater than mean(Z)
meanZ = np.mean(Z)
Zn = Zextreme[Zextreme > meanZ]
# now plot the Lorenz map Z_{n+1} vs Z{n}
plt.figure()
plt.plot(Zn[:-2],Zn[1:-1],'bo')
xlim = plt.gca().get_xlim()
plt.hold(True)
plt.plot(xlim,xlim,'k')
plt.xlabel('$Z_n$')
plt.ylabel('$Z_{n+1}$')
plt.title('Lorenz map: $r={}$, $\sigma={}$, $b={}$'.format(r,sigma,b))
plt.show()
Explanation: The Lorenz Map
The Lorenz map is simply a plot of the peak amplitude of the $n+1$ maximum $Z_{n+1}(t)$ as a function of the previous peak $Z_{n}$. Here we extract the maxima with a bit of discrete hackery where we
approximate the derivative of Z'
approximate extrema using intervals of the discrete solution with zero crossings of Z'
approximate local maxima using extrema with Z greater than mean(Z)
given the maxima $Z_n$,  plot $Z_{n+1}$ vs $Z_{n}$.
End of explanation
<END_TASK>