Hi, I am using package deSolve to run some ordinary differential equations
(ODE) as part of a mathematical modeling project. I have solved for the
following equilibrium states: Seq1<-a*(1-Neq1)/(f*Veq1+m+d)
Ceq1<-(f*Seq1*Veq1+g*Ieq1+r*(1-Neq1)-b1*Veq1*Ieq1)/(b2+m+d+g)
Ieq1<-(-b2*Ceq1)-r*(1-Neq1)/(b1*Veq1-g-u)
Veq1<-o*(Ceq1+Ieq1)/e I want to take the differential of both sides of the
equation and then solve for the inverse of the first as follows (the parameter
values are made up): library(deSolve)
rkMethod("rk45dp7")
CSIeq1<-function(t,yeq1,pars) {
with (as.list(c(yeq1,pars)),{
Neq1<-Seq1+Ceq1+Ieq1
dSeq1<-d[a*(1-Neq1)/(f*Veq1+m+d)]
dCeq1<-d[(f*Seq1*Veq1+g*Ieq1+r*(1-Neq1)-b1*Veq1*Ieq1)/(b2+m+d+g)]
dIeq1<-d[(-b2*Ceq1)-r*(1-Neq1)/(b1*Veq1-g-u)]
dVeq1<-d[o*(Ceq1+Ieq1)/e]
return(list(c(Seq1,Ceq1,Ieq1,Veq1),Neq1))
})
}pars <- c(a=0.1,
m=0.0005,
u=0.00005,
b1=0.7,
b2=0.2,
f=0.01,
g=0.4,
r=0.3,
o=20,
e=90,
d=0.01)
# initial conditions
yeq1 <- c(Seq1=0.99,Ceq1=0.01,Ieq1=0,Veq1=1)
t <-seq(0,365, by=50)
o1 <- ode(yeq1, t, CSIeq1, pars, method = rkMethod("rk45dp7"))1/o1$Seq1
*************************************************** The output gives the
following error message so I am wondering about the syntax of the differentials
such as dSeq1<-d[a*(1-Neq1)/(f*Veq1+m+d)]:
> library(deSolve)
> rkMethod("rk45dp7")
$ID
[1] "rk45dp7"$varstep
[1] TRUE$FSAL
[1] TRUE$A
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.00000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000
[2,] 0.20000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000
[3,] 0.07500000 0.225000 0.0000000 0.0000000 0.0000000 0.0000000
[4,] 0.97777778 -3.733333 3.5555556 0.0000000 0.0000000 0.0000000
[5,] 2.95259869 -11.595793 9.8228929 -0.2908093 0.0000000 0.0000000
[6,] 2.84627525 -10.757576 8.9064227 0.2784091 -0.2735313 0.0000000
[7,] 0.09114583 0.000000 0.4492363 0.6510417 -0.3223762 0.1309524$b1
[1] 0.08991319 0.00000000 0.45348907 0.61406250 -0.27151238 0.08904762
[7] 0.02500000$b2
[1] 0.09114583 0.00000000 0.44923630 0.65104167 -0.32237618 0.13095238
[7] 0.00000000$c
[1] 0.0000000 0.2000000 0.3000000 0.8000000 0.8888889 1.0000000 1.0000000$d
[1] -1.127018 0.000000 2.675424 -5.685527 3.521932 -1.767281
2.382469$densetype
[1] 1$stage
[1] 7$Qerr
[1] 4attr(,"class")
[1] "list" "rkMethod"
> CSIeq1<-function(t,yeq1,pars) {
+ with (as.list(c(yeq1,pars)),{
+ Neq1<-Seq1+Ceq1+Ieq1
+ dSeq1<-d[a*(1-Neq1)/(f*Veq1+m+d)]
+ dCeq1<-d[(f*Seq1*Veq1+g*Ieq1+r*(1-Neq1)-b1*Veq1*Ieq1)/(b2+m+d+g)]
+ dIeq1<-d[(-b2*Ceq1)-r*(1-Neq1)/(b1*Veq1-g-u)]
+ dVeq1<-d[o*(Ceq1+Ieq1)/e]
+ return(list(c(Seq1,Ceq1,Ieq1,Veq1),Neq1))
+ })
+ }
>
> pars <- c(a=0.1,
m=0.0005,
u=0.00005,
b1=0.7,
b2=0.2,
f=0.01,
g=0.4,
r=0.3,
o=20,
e=90,
d=0.01)
> # initial conditions
> yeq1 <- c(Seq1=0.99,Ceq1=0.01,Ieq1=0,Veq1=1)
> t <-seq(0,365, by=50)
> o1 <- ode(yeq1, t, CSIeq1, pars, method = rkMethod("rk45dp7"))
There were 50 or more warnings (use warnings() to see the first 50)
>
> 1/o1$Seq1
Error in o1$Seq1 : $ operator is invalid for atomic vectors
> warnings()
Warning messages:
1: In eval(expr, envir, enclos) : NAs introduced by coercion
******************************************** Any thoughts or suggestions would
be appreciated. If you run the program with only the differential on the left,
it runs just fine. Do I need to use a different 'R' package?
-- AAP
Anil A. Panackal, MD, ScM, FACP
Confidential E-Mail: This e-mail is intended only for the person or entity to
which it is addressed, and may contain information that is privileged,
confidential, or otherwise protected from disclosure. Dissemination,
distribution, or copying of this e-mail or the information herein by anyone
other than the intended recipient is prohibited. If you have received this
e-mail in error, please notify the sender by reply e-mail, and destroy the
original message and all copies.
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
