Anil A. Panackal MD
2012-Mar-02 19:33 UTC
[R] ?Syntax on Taking differential on both sides of the equation in 'R'
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$Seq1Error 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]]