odes.org
author Santosh G. Vattam <vattam.santosh@gmail.com>
Sun, 11 Apr 2010 02:30:55 +0530
changeset 42 303fe222243b
parent 14 4182c6c7e1c6
child 84 417992d2711e
permissions -rw-r--r--
Branches merged.
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
2
008c0edc6eac Added scripts for session-4 and session-6 of day-1.
Puneeth Chaganti <punchagan@gmail.com>
parents:
diff changeset
     1
* Solving ODEs
008c0edc6eac Added scripts for session-4 and session-6 of day-1.
Puneeth Chaganti <punchagan@gmail.com>
parents:
diff changeset
     2
*** Outline
008c0edc6eac Added scripts for session-4 and session-6 of day-1.
Puneeth Chaganti <punchagan@gmail.com>
parents:
diff changeset
     3
***** Introduction
008c0edc6eac Added scripts for session-4 and session-6 of day-1.
Puneeth Chaganti <punchagan@gmail.com>
parents:
diff changeset
     4
******* What are we going to do?
008c0edc6eac Added scripts for session-4 and session-6 of day-1.
Puneeth Chaganti <punchagan@gmail.com>
parents:
diff changeset
     5
******* How are we going to do?
008c0edc6eac Added scripts for session-4 and session-6 of day-1.
Puneeth Chaganti <punchagan@gmail.com>
parents:
diff changeset
     6
******* Arsenal Required
008c0edc6eac Added scripts for session-4 and session-6 of day-1.
Puneeth Chaganti <punchagan@gmail.com>
parents:
diff changeset
     7
********* working knowledge of arrays
008c0edc6eac Added scripts for session-4 and session-6 of day-1.
Puneeth Chaganti <punchagan@gmail.com>
parents:
diff changeset
     8
********* working knowledge of functions
008c0edc6eac Added scripts for session-4 and session-6 of day-1.
Puneeth Chaganti <punchagan@gmail.com>
parents:
diff changeset
     9
*** Script
14
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    10
    Welcome. 
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    11
    
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    12
    In this tutorial we shall look at solving Ordinary Differential Equations
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    13
    using odeints in Python.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    14
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    15
    Let's consider a classical problem of the spread of epidemic in a
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    16
    population.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    17
    This is given by dy/dt = ky(L-y) where L is the total population.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    18
    For our problem Let us use L=25000, k=0.00003.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    19
    Let the boundary condition be y(0)=250.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    20
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    21
    First of all run the magic command to import odeint to our program.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    22
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    23
    In []: from scipy.integrate import odeint
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    24
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    25
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    26
    For now just remember this as a command that does some magic to obtain
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    27
    the function odeint in to our program.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    28
    We will come back to the details of this command in subsequent sessions.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    29
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    30
    We can represent the given ODE as a Python function.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    31
    This function takes the dependent variable y and the independent variable t
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    32
    as arguments and returns the ODE.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    33
    Our function looks like this:
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    34
    (Showing the slide should be sufficient)
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    35
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    36
    In []: def epid(y, t):
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    37
      ....     k = 0.00003
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    38
      ....     L = 25000
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    39
      ....     return k*y*(L-y)
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    40
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    41
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    42
    Independent variable t can have be assigned the values in the interval of
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    43
    0 and 12 with 61 points using linspace:
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    44
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    45
    In []: t = linspace(0, 12, 61)
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    46
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    47
    Now obtaining the odeint of the ode we have already defined is as simple as
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    48
    calling the Python's odeint function which we imported:
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    49
    
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    50
    In []: y = odeint(epid, 250, t)
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    51
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    52
    We can plot the the values of y against t to get a graphical picture our ODE.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    53
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    54
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    55
    Let us move on to solving a system of two ordinary differential equations.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    56
    Here we shall take the example ODEs of a simple pendulum.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    57
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    58
    The equations can be written as a system of two first order ODEs
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    59
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    60
    d(theta)/dt = omega
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    61
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    62
    and
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    63
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    64
    d(omega)/dt = - g/L sin(theta)
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    65
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    66
    Let us define the boundary conditions as: at t = 0, 
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    67
    theta = theta 0 (10 degrees) and omega = 0
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    68
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    69
    Let us first define our system of equations as a Python function, pend_int.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    70
    As in the earlier case of single ODE we shall use odeint function of Python
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    71
    to solve this system of equations by passing pend_int to odeint.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    72
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    73
    pend_int is defined as shown:
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    74
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    75
    In []: def pend_int(initial, t):
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    76
      ....     theta = initial[0]
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    77
      ....     omega = initial[1]
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    78
      ....     g = 9.81
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    79
      ....     L = 0.2
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    80
      ....     f=[omega, -(g/L)*sin(theta)]
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    81
      ....     return f
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    82
      ....
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    83
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    84
    It takes two arguments. The first argument is a 2-tuple containing the two
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    85
    dependent variables in the system, theta and omega.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    86
    The second argument is the independent variable t.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    87
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    88
    In the function we assign theta and omega to first and second values of the
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    89
    initial argument respectively.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    90
    Acceleration due to gravity, as we know is 9.8 meter per second sqaure.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    91
    Let the length of the the pendulum be 0.2 meter.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    92
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    93
    We create a list, f, of two equations which corresponds to our two ODEs,
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    94
    that is d(theta)/dt = omega and d(omega)/dt = - g/L sin(theta).
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    95
    We return this list of equations f.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    96
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    97
    Now we can create a set of values for our time variable t over which we need
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    98
    to integrate our system of ODEs. Let us say,
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
    99
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   100
    In []: t = linspace(0, 20, 101)
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   101
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   102
    We shall assign the boundary conditions to the variable initial.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   103
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   104
    In []: initial = [10*2*pi/360, 0]
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   105
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   106
    Now solving this system is just a matter of calling the odeint function with
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   107
    the correct arguments.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   108
    So first let us import odeint function into our program using the magic
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   109
    import command
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   110
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   111
    In []: from scipy.integrate import odeint
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   112
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   113
    We can call ode_int as:
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   114
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   115
    In []: pend_sol = odeint(pend_int, initial,t)
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   116
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   117
    Plotting theta against t and omega against t we obtain the plots as shown
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   118
    in the slide.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   119
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   120
    Thus we come to the end of this session on solving ordinary differential
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   121
    equations in Python. Thanks for listening to this tutorial.
4182c6c7e1c6 Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents: 2
diff changeset
   122
2
008c0edc6eac Added scripts for session-4 and session-6 of day-1.
Puneeth Chaganti <punchagan@gmail.com>
parents:
diff changeset
   123
*** Notes