deSolve difference equations returning static output

Question

I am trying to model a series of difference equations in R using the deSolve package and method = "iteration", but the output is just returning the initial values with no changes, and I do not understand why. This is my code:

library(deSolve)

travelmodel <- function(t,y,parms) {
  with(as.list(c(y,parms)), {
    
    dImVillage <- Im*Village*(1-a) + ImField
    dImField <- Im*Village*a + Ih*(1-Im)*Village*a
    ImFlying <- ImField + Im*Village*a + Ih*(1-Im)*Village*a
    ImStaying <- 1- ImFlying
    
    dSmVillage <- (1-Im)*Village*(1-a) + SmField
    dSmField <- (1-Ih)*(1-Im)*Village*a 
    SmFlying <- SmField + (1-Ih)*(1-Im)*Village*a
    SmStaying <- 1- ImFlying
    
    Village = ImField + Im*Village*a + Ih*(1-Im)*Village*a + SmField + (1-Ih)*(1-Im)*Village*a
    
    list(c(ImVillage, ImField, ImFlying, ImStaying, SmVillage, SmField, SmFlying, SmStaying, Village))
  })
}  
parmstrav <- c(Ih = 0.5, 
               Im = 0.5, 
               a = 0.15)

times <- seq(0, 365, by = 1)   statetrav <- c(ImVillage = 0.9,
                 ImField = 0.1,
                 ImFlying = 0,
                 ImStaying = 0,
                 SmVillage = 0.9,
                 SmField = 0.1,
                 SmFlying = 0,
                 SmStaying = 0,
                 Village = 0)     
difference_output <- ode(y = statetrav, times = times, func = travelmodel, parms = parmstrav, method = "iteration")

I would really appreciate any help with this, thanks a lot.

Answer 1

Here a small example how to formulate a difference equation model by specifying either the differences or the new states in the model function. In the first case, method="euler" is used while method = "iteration" is needed in the second case. Both formulations below return identical results:

library(deSolve)

y0 <- c(n = 1, r = 1)
p <- c(k = 0.1)
times <- 0:100

## (1) formulation returns differences
differences <- function(t, y, p) {
  with(as.list(c(y, p)),{
    dn <- r * n
    dr <- -k * r

    list(c(dn, dr))
  })
}

out1 <- ode(y0, times, differences, p, method = "euler")
plot(out1)

## (2) formulation returns new states
new_states <- function(t, y, p) {
  with(as.list(c(y, p)),{
    n <- n + r * n
    r <- r - k * r
    list(c(n, r))
  })
}

out2 <- ode(y0, times, new_states, p, method = "iteration")
plot(out2)

## the two solutions are identical
summary(out1-out2)

Two additional notes:

• In case of difference equations (i.e. discrete time steps) the order of the equations matters.
• The first formulation can easily be solved as a continuous differential equation by using another solver, e.g. method = "lsoda".