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.



In [2]:

    
from IPython.html.widgets import interact, interactive, 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 [3]:

    
def lorentz_derivs(yvec, t, sigma, rho, beta):
    """Compute the the derivatives for the Lorentz system at yvec(t)."""
    # YOUR CODE HERE
    #     raise NotImplementedError()
    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 [4]:

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

    
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.
    
    """
    # YOUR CODE HERE
    #     raise NotImplementedError()
    t = np.linspace(0, max_time, 250*max_time)
    soln = odeint(lorentz_derivs, ic, t, args=(sigma, rho,beta), atol=1e-9, rtol=1e-9)
    return np.array(soln), np.array(t)



In [6]:

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

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

    
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.
    """
    # YOUR CODE HERE
#     raise NotImplementedError()
    plt.figure(figsize=(15,8)) 
    np.random.seed(1)
    r = []
    for i in range(0,10):
        data = (np.random.random(3)-0.5)*30.0
        r.append(solve_lorentz(data, max_time, sigma, rho, beta))
    for j in r:
        x = [p[0] for p in j[0]]
        z = [p[2] for p in j[0]]
        color = plt.cm.summer((x[0] + z[0]/60.0 - 0.5))
        plt.plot(x, z, color=color)
    plt.xlabel('$x(t)$')
    plt.ylabel('$z(t)$')
    plt.title('Lorentz system')



In [9]:

    
plot_lorentz()



In [10]:

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

    
# YOUR CODE HERE
# raise NotImplementedError()
w = interactive(plot_lorentz, max_time = (1,10,1), N = (1,50,1), sigma = (0.0,50.0,0.1), rho = (0.0, 50.0, 0.1), bata = fixed(8.0/3.0));
w

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

YOUR ANSWER HERE