1 * Solving ODEs

     2 *** Outline

     3 ***** Introduction

     4 ******* What are we going to do?

     5 ******* How are we going to do?

     6 ******* Arsenal Required

     7 ********* working knowledge of arrays

     8 ********* working knowledge of functions

     9 *** Script

    10     Welcome friends. 

    11     

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

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

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

    15     previous tutorials.

    16 

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

    18     population.

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

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

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

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

    23 

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

    25     

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

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

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

    29     So we run the magic command:

    30 

    31     In []: from scipy.integrate import odeint

    32 

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

    34 

    35     We can represent the given ODE as a Python function.

    36     This function takes the dependent variable y and the independent variable t

    37     as arguments and returns the ODE.

    38 

    39     Let us now define our function.

    40 

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

    42       ....     k = 0.00003

    43       ....     L = 250000

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

    45 

    46     Independent variable t can be assigned the values in the interval of

    47     0 and 12 with 61 points using linspace:

    48 

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

    50 

    51     Now obtaining the solution of the ode we defined, is as simple as

    52     calling the Python's odeint function which we just imported:

    53     

    54     In []: y = odeint(epid, 250, t)

    55 

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

    57 

    58     plot(y, t)

    59     Lets now close this plot and move on to solving ordinary differential equation of 

    60     second order.

    61     Here we shall take the example ODEs of a simple pendulum.

    62 

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

    64 

    65     d(theta)/dt = omega

    66     

    67     and

    68 

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

    70 

    71     Let us define the boundary conditions as: at t = 0, 

    72     theta = theta naught = 10 degrees and 

    73     omega = 0

    74 

    75     Let us first define our system of equations as a Python function, pend_int.

    76     As in the earlier case of single ODE we shall use odeint function of Python

    77     to solve this system of equations by passing pend_int to odeint.

    78 

    79     pend_int is defined as shown:

    80 

    81     In []: def pend_int(initial, t):

    82       ....     theta = initial[0]

    83       ....     omega = initial[1]

    84       ....     g = 9.81

    85       ....     L = 0.2

    86       ....     f=[omega, -(g/L)*sin(theta)]

    87       ....     return f

    88       ....

    89 

    90     It takes two arguments. The first argument itself containing two

    91     dependent variables in the system, theta and omega.

    92     The second argument is the independent variable t.

    93 

    94     In the function we assign the first and second values of the

    95     initial argument to theta and omega respectively.

    96     Acceleration due to gravity, as we know is 9.81 meter per second sqaure.

    97     Let the length of the the pendulum be 0.2 meter.

    98 

    99     We create a list, f, of two equations which corresponds to our two ODEs,

   100     that is d(theta)/dt = omega and d(omega)/dt = - g/L sin(theta).

   101     We return this list of equations f.

   102 

   103     Now we can create a set of values for our time variable t over which we need

   104     to integrate our system of ODEs. Let us say,

   105 

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

   107 

   108     We shall assign the boundary conditions to the variable initial.

   109 

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

   111 

   112     Now solving this system is just a matter of calling the odeint function with

   113     the correct arguments.

   114 

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

   116 

   117     In []: plot(pend_sol)

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

   119 

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

   121     differential equations of first and second order.

   122 

   123 *** Notes