Violet

10 Reputation

One Badge

1 years, 71 days

MaplePrimes Activity


These are replies submitted by Violet

@mmcdara thank you so much. I am now getting what I wanted to see on my graphs.

@mmcdara below is the rcode

Rcode.txtlibrary(deSolve)
spirs <- function(t, x, parms) {
with(as.list(c(parms, x)), {
dS=L-(1-u)*b*S*P+th*R-nu*S
dI=(1-u)*b*S*P-(nu+sg+v+m+a-ph)*I
dP= m*I-g1*P
dR=(v+a)*I-(nu+th)*R
output <- c(dS, dI,dP, dR)
list(output)
})
}
#the Initial values
#start<-c(S=0.8,P=0.08,I=0.12,R=0)
start<-c(S=1000,P=100,I=150,R=0)

## The parameters
parms <- c(L=0.05,b=0.002,ph=0.01,a=0.01,m=0.01,nu=0.001,th=0.005,u=0, v=0,sg=0.01, g1=0.03)

## vector of timesteps
times <- seq(0,60,1) ###odel run in 52 weeks*20 years
run_d<-ode(times=times, y=start, func=spirs,parms=parms)
summary(run_d)

matplot(run_d[,2:5], type="l")
plot(run_d[,2], lwd=3, col="black", ylim=c(0,1000), type="l", main="SIR Model", ylab="Population")
lines(run_d[,3], lwd=3,col="red")
lines(run_d[,4], lwd=3, col="blue")
lines(run_d[,5], lwd=3, col="green")
legend("topright",legend=c("S", "P","I", "R"), col=c("black","red", "blue", "green"), lty=c(1,1,1,1))

###R0
L=0.05;b=0.002;ph=0.01;a=0.01;m=0.01;nu=0.001;th=0.005;u=0; v=0;sg=0.01; g1=0.03
Ro=(1-u)*b*m*L/((g1*nu)*(m+v+a+sg+nu-ph))
Ro

Download Rcode.txt

 


@mmcdara below is my maple code

L := 0.5e-1:

Ro := (1-u)*b*m*L/(g1*nu*(nu+sg+v+m+a-ph))

.2479338843

(1)

with(plots):

plots[display]({H1, H2, H3, H4});

 

diff(S1(t), t) = L-(1-u)*b*S1(t)*P1(t)+th*R1(t)-nu*S1(t), diff(I1(t), t) = (1-u)*b*S1(t)*P1(t)-(nu+sg+v+m+a-ph)*I1(t), diff(P1(t), t) = m*I1(t)-g1*P1(t), diff(R1(t), t) = (a+v)*I1(t)-(nu+th)*R1(t)

diff(S1(t), t) = 0.1e-1-0.1e-2*S1(t)*P1(t)+0.1e-1*R1(t)-0.1e-2*S1(t), diff(I1(t), t) = 0.1e-2*S1(t)*P1(t)-.270*I1(t), diff(P1(t), t) = 0.1e-1*I1(t)-0.1e-1*P1(t), diff(R1(t), t) = .25*I1(t)-0.11e-1*R1(t)

(2)

``

with(plots):

L := 0.5e-1:

L := 0.5e-1:

L := 0.5e-1:

plots[display]({H3, H31, H32, H33});

 

restart

``


 

Download modie-work.mw

 

Page 1 of 1