The default numerical scheme used by the ode() function in the deSolve R package is the lsode method which implements the BDF and implicit Adams linear multisteps schemes to solve the ODE system. The integration is completed with a variable step size, that is controlled by estimating the local truncation error (LTE) at each step. In practice, the step size is adapted in such a way that the estimate of the LTE is smaller than some preset value rtol and atol. The default values in the ode method are:
ode(rtol = 1e-6, atol = 1e-6)
And can be adapted.
However, the actual evaluated LTE is always something above these values. Is there a way to extract this value from the solver?
An example system is given bellow:
rm(list = ls())
install.packages('deSolve')
library('deSolve')
# Example ODE system for the Lotka-V predator prey model
LVmod <- function(Time, State, Pars) {
with(as.list(c(State, Pars)), {
Ingestion <- rIng * Prey * Predator
GrowthPrey <- rGrow * Prey * (1 - Prey/K)
MortPredator <- rMort * Predator
dPrey <- GrowthPrey - Ingestion
dPredator <- Ingestion * assEff - MortPredator
return(list(c(dPrey, dPredator)))
})
}
# values of the parameters
pars <- c(rIng = 0.2, # /day, rate of ingestion
rGrow = 1.0, # /day, growth rate of prey
rMort = 0.2 , # /day, mortality rate of predator
assEff = 0.5, # -, assimilation efficiency
K = 10) # mmol/m3, carrying capacity
#initial values for the predator prey state variables
yini <- c(Prey = 1, Predator = 2)
#times at which the ode return state output
times <- seq(0, 200, by = 1)
# the solver using lsode
out <- ode(yini, times, LVmod, pars)
summary(out)
str(out)
deSolve [1:201, 1:3] 0 1 2 3 4 5 6 7 8 9 ...
- attr(*, "dimnames")=List of 2
..$ : NULL
..$ : chr [1:3] "time" "Prey" "Predator"
- attr(*, "istate")= int [1:21] 2 282 517 NA 1 1 0 52 22 NA ...
- attr(*, "rstate")= num [1:5] 1 1 201 0 143
- attr(*, "lengthvar")= int 2
- attr(*, "type")= chr "lsoda"
