

In [1]:

    
import numpy as np
import numpy.random as npr
import pylab as pl
from scipy.integrate import odeint
%matplotlib inline

scipy.integrate.odeint integrates ODEs of the form: $$ \newcommand{\y}{\mathbf{y}} \newcommand{\f}{\mathbf{f}} \frac{d \y}{d t} \equiv \dot{\y} = \f(\y,t)$$ where $\y \in \mathbb{R}^n$ with initial condition $$\y(t=t_0) = \y_0$$

Note that we can decompose $\y$ as the system: $$ \begin{align} \y & = (y_1,\dots,y_n)\\ \dot{y_1} & = f_1 (y_1,\dots,y_n,t)\\ \dot{y_2} & = f_2 (y_1,\dots,y_n,t)\\ &\dots\\ \dot{y_n} & = f_n (y_1,\dots,y_n,t) \end{align}$$

To use odeint, we define a function $\f : \mathbb{R}^{n+1} \to \mathbb{R}^n$, $\f: (\y,t) \mapsto \f(\y,t)$

1D example

Consider the one dimensional system: $$ \begin{align} n & = 1\\ y & \equiv \y\\ \dot{y} &= a y\\ a & = 2\\ y_0 &= 1\\ \end{align}$$ with exact solution $$ y(t) = y_0 \exp(at) = \exp(at)$$



In [2]:

    
# 1D example:
a=2
y0=1.0
f = lambda y,t:a*y



In [7]:

    
t = np.arange(0,6,0.02)
y = odeint(f,y0,t)[:,0]



In [10]:

    
pl.plot(t,y);

2D example

Consider the harmonic oscillator: $$ \ddot{x} + \omega^2 x = 0$$ where $\ddot{x}$ is the second-derivative of $x$ with respect to time, i.e. the time-derivative of $\dot{x}$ ($\ddot{x} \equiv \dot{y}$ where $y = \dot{x}$)

We need to re-write this as a system which contains no second-derivatives. We can do this by introducing a new variable $y \equiv \dot{x}$: $$ \begin{align} \dot{x} &= y\\ \dot{y} &= - \omega^2 x \end{align}$$

To re-write this as a system of first



In [10]:



In [12]:

    
def lorenz_sys(q,t):
    x = q[0]
    y = q[1]
    z = q[2]
    # sigma, rho and beta are global.
    f = [sigma * (y - x),
         rho*x - y - x*z,
         x*y - beta*z]
    return f



In [113]:

    
sigma = 10.0
rho = 28.0
beta = 10.0/3
y0 = np.ones(3)*10



In [114]:

    
t = np.arange(0,100,0.0001)
y = odeint(lorenz_sys,y0,t)



In [115]:

    
pl.plot(y[:,0])
pl.figure()
pl.plot(y[:,1])
pl.figure()
pl.plot(y[:,2])









    Out[115]:





[<matplotlib.lines.Line2D at 0x114a48fd0>]



In [116]:

    
pl.plot(y[:,0],y[:,1])









    Out[116]:





[<matplotlib.lines.Line2D at 0x1148eca50>]



In [117]:

    
from mpl_toolkits.mplot3d import Axes3D



In [118]:

    
fig = pl.figure()
ax = Axes3D(fig)
ax.plot(y[:,0], y[:,1], y[:,2])
pl.xlabel('x')
pl.ylabel('y')
pl.show()



In [119]:

    
fig = pl.figure()
ax = Axes3D(fig)
ax.plot(y[20000:,0], y[20000:,1], y[20000:,2])
pl.xlabel('x')
pl.ylabel('y')
pl.show()



In [95]:



In [ ]: