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 t,soln

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]:
colors = plt.cm.hot(1)

In [7]:
colors


Out[7]:
(0.051894844054434848, 0.0, 0.0, 1.0)

In [8]:
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 [33]:
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.
    """
    colors = plt.cm.hot(np.linspace(0,1,N))
    np.random.seed(1)
    for i in range(N):
        ic = (np.random.random(3)-.5)*30
        t,soln = solve_lorentz(ic, max_time, sigma, rho, beta)
        x = [e[0] for e in soln]
        y = [e[1] for e in soln]
        z = [e[2] for e in soln]
        plt.plot(x,z,color=colors[i])
        plt.title('Lorenz System for Multiple Trajectories')
        plt.xlabel('X Position')
        plt.ylabel('Y Position')

In [34]:
plot_lorentz()



In [35]:
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 [36]:
interact(plot_lorentz,max_time=(1,11,1),N=(1,51,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:

As I vary $\sigma$ to high values, the curves become more spread out and symmetric, but if sigma is too low the curves become lines and the system stops representing chaotic motion.

As I increase $\rho$ the number of circular curves increases. So each trajectory spins around more times.

As I decrease $\beta$ the spirals die very quickly

I worked with Hunter, Jessica, and Brett.