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

    
np.random.seed(1)
g=[]
h=[]
f=[]
for i in range(5):
    rnd = np.random.random(size=3)
    a,b,c = 30*rnd - 15
    g.append(a)
    h.append(b)
    f.append(c)
g,h,f









    Out[7]:





([-2.4893398589227793,
  -5.9300228210448065,
  -9.4121936586698727,
  1.1645020201007092,
  -8.8664325080544764],
 [6.6097348032647432,
  -10.597323275486609,
  -4.6331781887085679,
  -2.4241645679011565,
  11.343523091728361],
 [-14.996568755479654,
  -12.229842156936066,
  -3.096975773079901,
  5.5565850119027864,
  -14.178372204062216])



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.
    """
    np.random.seed(1)
    colors = plt.cm.hot(np.linspace(0,1,N))
    f = plt.figure(figsize=(7,7))
    for i in range(N):
        ic = 30*np.random.random(size=3) - 15
        soln, t = solve_lorentz(ic, max_time, sigma, rho, beta)
        plt.plot(soln[:,0], soln[:,2], color=colors[i])
    plt.xlabel('x(t)')
    plt.ylabel('z(t)')
    plt.title('Lorenz System: x(t) vs. z(t)')
    plt.ylim(-20,110)
    plt.xlim(-60,60)



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

    
interact(plot_lorentz, N=(1,50), max_time=(1,10), 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$ = 25 gives a butterfly look to the graph. A low $\sigma$ value gives more V-shaped graph, with $\sigma$ = 0 creating vertical lines of varying lengths and positions (dependent on ic). $\rho$ affects the overall size of the graph, with larger $\rho$ corresponding to a larger graph.