View source for Python-Old-Version/C3/Ordinary-Differential-Equations /English

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.

 In []: 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.