Ordinary Differential Equations Exercise 1

Imports



In [26]:

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

    
def lorentz_derivs(yvec, t, sigma, rho, beta):
    """Compute the the derivatives for the Lorentz system at yvec(t)."""
    x,y,z = yvec[0], yvec[1], yvec[2]
    sol = [0,0,0]
    sol[0] = sigma*(y-x)
    sol[1] = x*(rho-z)-y
    sol[2] = x*y - beta*z
    return sol



In [28]:

    
np.zeros((2,3))









    Out[28]:





array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])



In [29]:

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

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



In [31]:

    
solve_lorentz(0,0,0)









    Out[31]:





(array([], shape=(0, 1), dtype=float64), array([], dtype=float64))



In [32]:

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

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

    
np.random.seed(1)
30 * np.random.random_sample((3,)) - 15
30 * np.random.random_sample((3,)) - 15









    Out[34]:





array([ -5.93002282, -10.59732328, -12.22984216])



In [35]:

    
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))
    for i in range(N):
        ic = 30 * np.random.random_sample((3,)) - 15
        ans, t = solve_lorentz(ic, max_time, sigma, rho, beta)
        c = colors[i]
        plt.plot(ans[0],ans[2],color =c )
        
    plt.xlim(-15,15)
    plt.ylim(-15,20)
    plt.xlabel('x(t)')
    plt.ylabel('z(t)')
    plt.title('Lorentz System')



In [36]:

    
plot_lorentz()



In [37]:

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

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

An increase in the sigma causes more lines to have on average a lower slope. A increase in rho causes more chaotic behavior. We kept beta constant so i am not sure of the change.