Revision history [back]
 | 1 | initial version |
Here is one solution using ode_solver.
var('S,E,I,W')
S0=2600
E0=1
I0=1
W0=200
bR= 37
k= 8760
m= .0142
r= .9
B= .0255*(1/W0)
A= 1
Y=100000
mW=26
timedata=srange(0,10,0.05)
f=[bR-B*S*W-m*S+k*I,B*S*W-A*E-m*E,A*E-k*I-m*I,Y*I-mW*W]
soln=desolve_odeint(f,[S0,E0,I0,W0],timedata,[S,E,I,W])
splot=line(zip(timedata,[s[0] for s in soln]),color='green')
eplot=line(zip(timedata,[s[1] for s in soln]),color='blue')
iplot=line(zip(timedata,[s[2] for s in soln]),color='red')
wplot=line(zip(timedata,[s[3] for s in soln]),color='black')
show(splot)
show(eplot)
show(iplot
show(wplot)
| | 2 | No.2 Revision |
Here is one solution using ode_solverdesolve_odeint
var('S,E,I,W')
S0=2600
E0=1
I0=1
W0=200
bR= 37
k= 8760
m= .0142
r= .9
B= .0255*(1/W0)
A= 1
Y=100000
mW=26
timedata=srange(0,10,0.05)
f=[bR-B*S*W-m*S+k*I,B*S*W-A*E-m*E,A*E-k*I-m*I,Y*I-mW*W]
soln=desolve_odeint(f,[S0,E0,I0,W0],timedata,[S,E,I,W])
splot=line(zip(timedata,[s[0] for s in soln]),color='green')
eplot=line(zip(timedata,[s[1] for s in soln]),color='blue')
iplot=line(zip(timedata,[s[2] for s in soln]),color='red')
wplot=line(zip(timedata,[s[3] for s in soln]),color='black')
show(splot)
show(eplot)
show(iplot
show(wplot)
