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