[Tutor] System of ODEs Question

[Eike Welk]

Hello Eric!


On Thursday 03.02.2011 08:05:13 Eric Lofgren wrote:

> def eq_system(PopIn,x):
>     '''Defining SIR System of Equations'''
>     #Creating an array of equations
>     Eqs= np.zeros((3))
>     Eqs[0]= -beta * PopIn[0]*PopIn[1]
>     Eqs[1]= beta * PopIn[0]*PopIn[1] - gamma*PopIn[1]
>     Eqs[2]= gamma*PopIn[1]
>     return Eqs
> 
> SIR = spi.odeint(eq_system, PopIn, t_interval)
> print SIR

> The part that is confusing me is defining the function. Namely, it seems to
> need two arguments - the first needs to be the initial states of the
> population being modeled, but the second...it doesn't seem to matter what
> it is. Originally, I didn't include it at all, and the program was very,
> very unhappy - looking back at the example, it had a second argument, so I
> put one in, and magically, it worked. But it can be x, it can be t, it can
> be anything.

The meaning of the arguments of the system function (`eq_system` in your code) 
is as follows:

1. (`PopIn` in your code) the current state of the system
2. (`x`     in your code) the current time

The arguments can have any namem, you are free to choose what you think is 
most appropriate. 

The documentation is here:
http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

In your system function the time is not used because the system is time 
invariant. But in a more advanced model you might use the time too, for 
example the virus particles might survive longer outside the human body in 
winter (`beta` a function of time).

I think a good way to write system functions is like this:

def sir_system(states, time):
    '''Defining SIR System of Equations'''
    s, i, r = states
    ds_dt = -beta * s * i
    di_dt = beta * s * i - gamma * i
    dr_dt = gamma * i
    return np.array([ds_dt, di_dt, dr_dt])


Eike.