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

    
def lorentz_derivs(yvec, t, sigma, rho, beta):
    """Compute the the derivatives for the Lorentz system at yvec(t)."""
    # YOUR CODE HERE
    raise NotImplementedError()



In [ ]:

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

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



In [ ]:

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

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



In [ ]:

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



In [ ]:

    
plot_lorentz()



In [ ]:

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

    
# YOUR CODE HERE
raise NotImplementedError()

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

YOUR ANSWER HERE