0

Very new user here. I am trying to use lsoda to solve differential equations stratified into two layers (as denoted by the for(s in 1:2) loop).

When running this full code, I keep getting the error message

object 'N' not found

no matter where or how I try to define N.

Can anyone help spot the error or advise on what I'm doing wrong? Thanks in advance.

R code:

library(deSolve)

Dyn <- function(t, var,par) {

  with(as.list(c(par, var)), {

    for(s in 1:2){ 

      #Derivatives
      dX[s] <- mu*N[s] - sigma*X[s] - (c[s]*beta*(InD[s] +ID[s]+ IdT[s])/N[s])*X[s] - mu*X[s]

      dXint[s] <- sigma*X[s] - (1-omega)*(c[s]*beta*(InD[s] +ID[s]+ IdT[s])/N[s])*Xint[s] - mu*Xprep[s] 

      dInD[s] <- (c[s]*beta*(InD[s] +ID[s]+ IdT[s])/N[s])*X[s] - psi*InD[s]- mu*InD[s]

      dID[s] <- (1-omega)*(c[s]*beta*(InD[s] +ID[s]+ IdT[s]) /N[s])*Xint[s]+ psi*InD[s]- mu*ID[s]

      N[s] <- X[s]+Xint[s]+InD[s]+ID[s]

      diffs <- c(dX[s], dXint[s], dInD[s], dID[s], N[s])}

    return(list(diffs))

  })}

#Defining parameter and initial values 

par <- c(mu=0.033, sigma=0.29, beta=0.40, c=c(2, 30), Ctot=1773600, N=c(332550, 36950), psi=0.022, omega=0.44)

init <- c(X=c(332550,36950), Xint=c(0,0), InD=c(1,1), ID=c(0,0))

t <- seq(0, 30, by=0.1) 

#Numerical solution#

Hom.sol <- lsoda(init, t, Dyn,par)
user6767735
  • 3
  • 3
  • I dont see where you've defined ``N`` -provide us a full reproducible example. For help on that see [here](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) – Cyrus Mohammadian Aug 28 '16 at 20:51

1 Answers1

0

I think you are mixing up parameters and variables. N seems to be defined as a parameter par with dimension 2. However, in your model definition you are updating N with dimension 1.

csgillespie
  • 59,189
  • 14
  • 150
  • 185
  • Thanks for spotting this. I fixed the code (see above) so that N now has 2 dimensions also ---by adding N[s]. I still get the same error message--"object 'N' not found". N[s] is the total of all other compartments in that strata, as now defined. Any thoughts to what I'm missing? Thx. – user6767735 Aug 31 '16 at 15:39