author | Shantanu <shantanu@fossee.in> |
Mon, 19 Apr 2010 12:45:13 +0530 | |
changeset 84 | 417992d2711e |
parent 14 | 4182c6c7e1c6 |
child 86 | a63a14de8584 |
permissions | -rw-r--r-- |
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 |
|
84 | 21 |
Lets start ipython -pylab interpreter. |
22 |
||
23 |
As we saw in one of earlier session, sometime pylab wont 'import' all |
|
24 |
packages. For solving 'ordinary differential equations' also we shall |
|
25 |
import 'odeint' function which is part SciPy package. So we run the |
|
26 |
magic command: |
|
14
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
27 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
28 |
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
|
29 |
|
84 | 30 |
# For now just remember this as a command that does some magic to obtain |
31 |
# the function odeint in to our program. |
|
32 |
We will cover more details regarding 'import' in subsequent sessions. |
|
14
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
33 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
34 |
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
|
35 |
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
|
36 |
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
|
37 |
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
|
38 |
(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
|
39 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
40 |
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
|
41 |
.... k = 0.00003 |
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
42 |
.... L = 25000 |
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
43 |
.... 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
|
44 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
45 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
46 |
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
|
47 |
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
|
48 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
49 |
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
|
50 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
51 |
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
|
52 |
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
|
53 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
54 |
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
|
55 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
56 |
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
|
57 |
|
84 | 58 |
plot(y, t) |
59 |
Lets close this plot and move on to solving ordinary differential equation of |
|
60 |
second order. |
|
14
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
61 |
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
|
62 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
63 |
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
|
64 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
65 |
d(theta)/dt = omega |
84 | 66 |
|
14
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
67 |
and |
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 |
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
|
70 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
71 |
Let us define the boundary conditions as: at t = 0, |
84 | 72 |
theta = theta naught = 10 degrees and |
73 |
omega = 0 |
|
14
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 |
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
|
76 |
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
|
77 |
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
|
78 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
79 |
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
|
80 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
81 |
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
|
82 |
.... theta = initial[0] |
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
83 |
.... omega = initial[1] |
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
84 |
.... g = 9.81 |
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
85 |
.... L = 0.2 |
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
86 |
.... 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
|
87 |
.... return f |
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
88 |
.... |
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
89 |
|
84 | 90 |
It takes two arguments. The first argument itself containing two |
14
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
91 |
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
|
92 |
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
|
93 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
94 |
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
|
95 |
initial argument respectively. |
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
96 |
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
|
97 |
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
|
98 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
99 |
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
|
100 |
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
|
101 |
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
|
102 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
103 |
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
|
104 |
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
|
105 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
106 |
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
|
107 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
108 |
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
|
109 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
110 |
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
|
111 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
112 |
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
|
113 |
the correct arguments. |
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 |
|
84 | 117 |
In []: plot(pend_sol[0], t) plot theta against t |
118 |
In []: plot(pend_sol[1], t) will plot omega against t |
|
14
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
119 |
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
|
120 |
in the slide. |
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
121 |
|
4182c6c7e1c6
Added script on ODEs, session 6 of our workshop.
Madhusudan.C.S <madhusudancs@gmail.com>
parents:
2
diff
changeset
|
122 |
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
|
123 |
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
|
124 |
|
2
008c0edc6eac
Added scripts for session-4 and session-6 of day-1.
Puneeth Chaganti <punchagan@gmail.com>
parents:
diff
changeset
|
125 |
*** Notes |