notebook.community

Edit and run

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) 




In [ ]:



    


%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

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





In [ ]:



    


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()

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





In [ ]:



    


# 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)

Damped Oscillation $r=10$

Now we increase the Rayleigh number to $r=10$ which admits two steady solutions depending on initial condition.





In [ ]:



    


# 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)

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.





In [ ]:



    


# 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)

Limit Cycle at large Rayleigh number

Now we increase the Rayleigh number to $r=350$ and the solution goes to a periodic limit cycle





In [ ]:



    


# 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)

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





In [ ]:



    


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()



    











In [ ]:

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.





In [ ]:



    


# 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()

Estimating the Liapunov exponent

Now we will calculate and plot the distance between the two trajectories on a log plot





In [ ]:



    


# 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()

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.





In [ ]:



    


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()

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

  1. approximate the derivative of Z'
  2. approximate extrema using intervals of the discrete solution with zero crossings of Z'
  3. approximate local maxima using extrema with Z greater than mean(Z)
  4. given the maxima $Z_n$, plot $Z_{n+1}$ vs $Z_{n}$.




In [ ]:



    


# 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()



    











In [ ]:

Content source: mspieg/dynamical-systems

Similar notebooks:

notebook.community | gallery | about