Ordinary Differential Equations Exercise 1

Imports



In [1]:

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









    



:0: FutureWarning: IPython widgets are experimental and may change in the future.

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

    
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 np.array([dx, dy, dz])



In [3]:

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

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



In [5]:

    
assert True # leave this to grade solve_lorenz

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

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

    
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)
    colors = plt.cm.hot(np.linspace(0,1,N))
    x = []
    y = []
    z = []
    for i in range(N):
        x.append(np.random.uniform(-15.0, 15.0))
        y.append(np.random.uniform(-15.0, 15.0))
        z.append(np.random.uniform(-15.0, 15.0))
    X = np.array(x)
    Y = np.array(y)
    Z = np.array(z)
    ic = np.transpose(np.vstack((X,Y,Z)))
    for i in range(len(X)):
        soln = solve_lorentz(ic[i], max_time, sigma, rho, beta)[0]
        soln2 = np.transpose(soln)
        plt.plot(soln2[0], soln2[2], color=colors[i])
    plt.xlabel("x(t)")
    plt.ylabel("z(t)")
    plt.title("Lorentz Plot")
    ax = plt.gca()
    ax.spines['top'].set_color('white')
    ax.spines['right'].set_color('white')
    plt.show()



In [8]:

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

    
interact(plot_lorentz, max_time = (1,10), N=(1,50), 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:

$\sigma$ increases the number of spirals. $\rho$ changes the distance of the spirals from each other. $\beta$ alters the values of spiral intercepts.



In [ ]: