* 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 friends. 

    In this tutorial we shall look at solving Ordinary Differential Equations,

    ODE henceforth using odeint in Python. In this tutorial we shall be using

    the concepts of arrays, functions and lists which we have covered in 

    previous tutorials.

    Let's consider the classic problem of the spread of an epidemic in a

    population.

    This is given by the ordinary differential equation dy/dt = ky(L-y) 

    where L is the total population and k is an arbitrary constant. For our

    problem Let us use L=250000, k=0.00003.

    Let the boundary condition be y(0)=250.

    Let's now fire up IPython by typing ipython -pylab interpreter.    

    As we saw in one of the earlier sessions, sometimes we will need more than 

    pylab to get our job done. For solving 'ordinary differential equations'

    we shall have to import 'odeint' function, which is a part of the SciPy package.

    So we run the magic command:

    In []: from scipy.integrate import odeint

    The details regarding `import' shall be covered in a subsequent tutorial.

    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.

    Let us now define our function.

    In []: def epid(y, t):

      ....     k = 0.00003

      ....     L = 250000

      ....     return k*y*(L-y)

    Independent variable t can 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 solution of the ode we defined, is as simple as

    calling the Python's odeint function which we just 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 now 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 the first and second values of the

    initial argument to theta and omega respectively.

    Acceleration due to gravity, as we know is 9.81 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)

    This gives us 2 plots. The green plot is omega vs t and the blue is theta vs t.

    Thus we come to the end of this tutorial. In this tutorial we have learnt how to solve ordinary

    differential equations of first and second order.

*** Notes