Ordinary Differential Equations Exercise 1

Imports



In [156]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from IPython.html.widgets import interact, fixed

Lorenz system

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 [157]:

    
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((dx,dy,dz))



In [158]:

    
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 [159]:

    
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,max_time*250)
    soln=odeint(lorentz_derivs,ic,t, args=(sigma,rho,beta))
    
    return np.array(soln),np.array(t)



In [ ]:



In [160]:

    
assert True # leave this to grade solve_lorenz



In [ ]:

Write a function plot_lorentz that:

Solves the Lorenz system for 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.
Plot $[x(t),z(t)]$ using a line to show each trajectory.
Color each line using the hot colormap from Matplotlib.
Label your plot and choose an appropriate x and y limit.

The following cell shows how to generate colors that can be used for the lines:



In [161]:

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









    



[ 0.0416  0.      0.      1.    ]
[ 0.70047002  0.          0.          1.        ]
[ 1.         0.3593141  0.         1.       ]
[ 1.          1.          0.02720491  1.        ]
[ 1.  1.  1.  1.]



In [168]:

    
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)
    
    
    plt.figure(figsize=(10,10))
    for i in range(N):                               #Uniform random samples
        x=np.random.uniform(-15,15)
        y=np.random.uniform(-15,15)
        z=np.random.uniform(-15,15)
        colors = plt.cm.hot(np.linspace(0,1,N))
        solutionarray, timearray=solve_lorentz((x,y,z),max_time,sigma,rho,beta)
        plt.plot(solutionarray[:,0], solutionarray[:,2], color=colors[i]) # plotting x column vs z column
        
    print(solutionarray)
    
    plt.xlabel('x(t)')
    plt.ylabel('z(t)')
    plt.set_cmap('hot')  # whether or not I have plt.set_cmap doesn't seem to matter for some reason..                       
    plt.title('Lorentz Solutions')
    plt.box(False)



In [169]:

    
plot_lorentz()









    



[[-13.8283565   -9.90508741  11.3442751 ]
 [-13.69156062 -10.76901538  11.79031188]
 [-13.59331797 -11.59647655  12.27272759]
 ..., 
 [ 16.1306741   19.35847338  34.30881879]
 [ 16.2473985   18.84463462  35.17646154]
 [ 16.33825745  18.27471907  36.00761134]]



In [164]:

    
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 [165]:

    
interact(plot_lorentz, max_time=([1,10]), N=([1,50]), sigma=([0.0,50.0]), rho=([0.0,50.0]), beta=fixed(8/3))









    



[[-13.8283565   -9.90508741  11.3442751 ]
 [-13.69156062 -10.76901538  11.79031188]
 [-13.59331797 -11.59647655  12.27272759]
 ..., 
 [ 16.1306741   19.35847338  34.30881879]
 [ 16.2473985   18.84463462  35.17646154]
 [ 16.33825745  18.27471907  36.00761134]]






    Out[165]:





<function __main__.plot_lorentz>

Describe the different behaviors you observe as you vary the parameters $\sigma$, $\rho$ and $\beta$ of the system:

-increasing/decreasing $\sigma$ appears to increase/decrease the overall distances from each point to another. So increasing $\sigma$ expands our system while decreasing $\sigma$ contracts our system.

-Increases/decreasing $\rho$ appears to increases/decreases the number of spiraling the system does.

-$\beta$ is fixed.



In [ ]: