Sorry, this content is no longer available

Ask Your Question
0

graphing ODE's using euler's method

asked 12 years ago

anonymous user

Anonymous

updated 10 years ago

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
Preview: (hide)

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 ( 12 years ago )

2 Answers

Sort by » oldest newest most voted
1

answered 12 years ago

calc314 gravatar image

updated 12 years ago

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)
Preview: (hide)
link
0

answered 12 years ago

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'))
Preview: (hide)
link

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 ( 12 years ago )

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: 12 years ago

Seen: 922 times

Last updated: Dec 05 '12