So I'm in the process of learning Python, and have been working on a program to
solve a very simple S-I-R model, used to study infectious diseases. Basically,
it's just a smallish set of differential equations which need to be numerically
integrated over time.
Working off of a template program, I came up with the following:
---
#Python implementation of continuous SIR model
#Import Necessary Modules
import numpy as np
import pylab as pl
import scipy.integrate as spi
#Parameter Values
S0 = 0.99999
I0 = 0.00001
R0 = 0.0
PopIn= (S0, I0, R0)
beta=1.4547
gamma=1/7.
t_end = 70
t_start = 1
t_step = 1
t_interval = np.arange(t_start, t_end, t_step)
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
#Plot Everything
#This is really ugly, but works for now
pl.plot(SIR[:,0])
pl.plot(SIR[:,1])
pl.plot(SIR[:,2])
pl.show()
---
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.
Which raises the question: What is it *doing*? Honestly, I'm baffled - can
anyone point me in the right direction?
Thanks,
Eric
_______________________________________________
Tutor maillist - [email protected]
To unsubscribe or change subscription options:
http://mail.python.org/mailman/listinfo/tutor
