# graphing ODE's using euler's method

Anonymous

Im trying to plot a graph based on the change in susceptible, exposed and infected people due to guinea worm, included is a plot of the worm population W. This is an attempt to at least replicate some of what is from this mathematical model

hopefully using the Impulsive differential equations but right now I've been relentlessly trying to graph these equations on a reasonable axis over ten years but for some reason it doesn't. Plus the y axis is completely off from what it should be, does anyone know what is wrong with the code? Thanks for any responses

timedata=[]
Sdata=[]
Edata=[]
Idata=[]
Wdata=[]

t=0
S=2600
E=1
I=1
W=200
bR= 37
k= 8760
m= .0142
r= .9
B= .0255*(1/W)
A= 1
Y=100000
mW=26
dt=.1

timedata.append(t)
Sdata.append(S)
Edata.append(E)
Idata.append(I)
Wdata.append(W)

T=10
Tfinal=(T/dt)
for i in range(0,Tfinal):
t= t+dt
Sprime=bR-B*S*W-m*S+k*I
Eprime=B*S*W-A*E-m*E
Iprime=A*E-k*I-m*I
Wprime=Y*I-mW*W
S= S+(Sprime*dt)
E= E+(Eprime*dt)
I= I+(Iprime*dt)
W= W+(Wprime*dt)
timedata.append(t)
Sdata.append(S)
Edata.append(E)
Idata.append(I)
Wdata.append(W)

Splot=list_plot(zip(timedata,Sdata),color='green',plotjoined=True)
Eplot=list_plot(zip(timedata,Edata),color='orange',plotjoined=True)
Iplot=list_plot(zip(timedata,Idata),color='black',plotjoined=True)

Splot

edit retag close merge delete

I suggest you take a closer look at your data or how you implemented your algorithm. I'm seeing stuff like 9.38387399424439e2962 in just the first bit of the Sdata, which I suspect is several (!) orders of magnitude bigger than what you're looking for...

( 2012-12-05 22:11:47 +0200 )edit

Sort by » oldest newest most voted

Here is one solution using desolve_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)

more

Here's an option using ode_solve.

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

T = ode_solver()
T.function=lambda t,y:[bR-B*y[0]*y[3]-m*y[0]+k*y[2],B*y[0]*y[3]-A*y[1]-m*y[1],A*y[1]-k*y[2]-m*y[2],Y*y[2]-mW*y[3]]
T.ode_solve(y_0=[S0,E0,I0,W0],t_span=[0,10],num_points=500)

a=T.solution
tssoln=[[soln[0],soln[1][0]] for soln in a]
tesoln=[[soln[0],soln[1][1]] for soln in a]
tisoln=[[soln[0],soln[1][2]] for soln in a]
twsoln=[[soln[0],soln[1][3]] for soln in a]

show(line(tssoln,color='green'))
show(line(tesoln,color='blue'))
show(line(tisoln,color='red'))
show(line(twsoln,color='black'))

more

You know, maybe these should go in the standard documentation as more complex examples, or be added to one of the PREP tutorials...

( 2012-12-06 08:12:57 +0200 )edit