* Solving ODEs
*** Outline
***** Introduction
******* What are we going to do?
******* How are we going to do?
******* Arsenal Required
********* working knowledge of arrays
********* working knowledge of functions
*** Script
    Welcome. 
    
    In this tutorial we shall look at solving Ordinary Differential Equations
    using odeints in Python.

    Let's consider a classical problem of the spread of epidemic in a
    population.
    This is given by dy/dt = ky(L-y) where L is the total population.
    For our problem Let us use L=25000, k=0.00003.
    Let the boundary condition be y(0)=250.

    Lets start ipython -pylab interpreter.    
    
    As we saw in one of earlier session, sometime pylab wont 'import' all
    packages. For solving 'ordinary differential equations' also we shall
    import 'odeint' function which is part SciPy package. So we run the 
    magic command:

    In []: from scipy.integrate import odeint

    # For now just remember this as a command that does some magic to obtain
    # the function odeint in to our program.
    We will cover more details regarding 'import' in subsequent sessions.

    We can represent the given ODE as a Python function.
    This function takes the dependent variable y and the independent variable t
    as arguments and returns the ODE.
    Our function looks like this:
    (Showing the slide should be sufficient)

    In []: def epid(y, t):
      ....     k = 0.00003
      ....     L = 25000
      ....     return k*y*(L-y)


    Independent variable t can have be assigned the values in the interval of
    0 and 12 with 61 points using linspace:

    In []: t = linspace(0, 12, 61)

    Now obtaining the odeint of the ode we have already defined is as simple as
    calling the Python's odeint function which we imported:
    
    In []: y = odeint(epid, 250, t)

    We can plot the the values of y against t to get a graphical picture our ODE.

    plot(y, t)
    Lets close this plot and move on to solving ordinary differential equation of 
    second order.
    Here we shall take the example ODEs of a simple pendulum.

    The equations can be written as a system of two first order ODEs

    d(theta)/dt = omega
    
    and

    d(omega)/dt = - g/L sin(theta)

    Let us define the boundary conditions as: at t = 0, 
    theta = theta naught = 10 degrees and 
    omega = 0

    Let us first define our system of equations as a Python function, pend_int.
    As in the earlier case of single ODE we shall use odeint function of Python
    to solve this system of equations by passing pend_int to odeint.

    pend_int is defined as shown:

    In []: def pend_int(initial, t):
      ....     theta = initial[0]
      ....     omega = initial[1]
      ....     g = 9.81
      ....     L = 0.2
      ....     f=[omega, -(g/L)*sin(theta)]
      ....     return f
      ....

    It takes two arguments. The first argument itself containing two
    dependent variables in the system, theta and omega.
    The second argument is the independent variable t.

    In the function we assign theta and omega to first and second values of the
    initial argument respectively.
    Acceleration due to gravity, as we know is 9.8 meter per second sqaure.
    Let the length of the the pendulum be 0.2 meter.

    We create a list, f, of two equations which corresponds to our two ODEs,
    that is d(theta)/dt = omega and d(omega)/dt = - g/L sin(theta).
    We return this list of equations f.

    Now we can create a set of values for our time variable t over which we need
    to integrate our system of ODEs. Let us say,

    In []: t = linspace(0, 20, 101)

    We shall assign the boundary conditions to the variable initial.

    In []: initial = [10*2*pi/360, 0]

    Now solving this system is just a matter of calling the odeint function with
    the correct arguments.

    In []: pend_sol = odeint(pend_int, initial,t)

    In []: plot(pend_sol[0], t) plot theta against t
    In []: plot(pend_sol[1], t) will plot omega against t
    Plotting theta against t and omega against t we obtain the plots as shown
    in the slide.

    Thus we come to the end of this session on solving ordinary differential
    equations in Python. Thanks for listening to this tutorial.

*** Notes