Ask Your Question
0

graphing ODE's using euler's method

asked 2012-12-05 14:59:18 -0500

anonymous user

Anonymous

updated 2015-01-13 11:20:17 -0500

FrédéricC gravatar image

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

(http://mysite.science.uottawa.ca/rsmi... ),

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 flag offensive close merge delete

Comments

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...

kcrisman gravatar imagekcrisman ( 2012-12-05 15:11:47 -0500 )edit

2 answers

Sort by » oldest newest most voted
1

answered 2012-12-05 15:50:57 -0500

calc314 gravatar image

updated 2012-12-05 16:08:34 -0500

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)
edit flag offensive delete link more
0

answered 2012-12-05 16:07:55 -0500

calc314 gravatar image

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'))
edit flag offensive delete link more

Comments

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

kcrisman gravatar imagekcrisman ( 2012-12-06 01:12:57 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2012-12-05 14:59:18 -0500

Seen: 207 times

Last updated: Dec 05 '12