Извлечение значения локальной ошибки усечения (LTE) на каждом шаге интеграции при использовании функции ode () из R-пакета deSolve - PullRequest
Новая
       43

Извлечение значения локальной ошибки усечения (LTE) на каждом шаге интеграции при использовании функции ode () из R-пакета deSolve

1 голос
/ 09 апреля 2019

Цифровой схемой по умолчанию, используемой функцией ode () в пакете deSolve R, является метод lsode, который реализует линейные многошаговые схемы BDF и неявных Адамса для решения системы ODE. Интеграция завершается с переменным размером шага, который контролируется путем оценки локальной ошибки усечения (LTE) на каждом шаге. На практике размер шага адаптируется таким образом, что оценка LTE меньше, чем некоторое заданное значение rtol и atol. Значения по умолчанию в методе ode: 

ode(rtol = 1e-6, atol = 1e-6)

И может быть адаптирован.

Однако фактическая оцененная LTE всегда находится выше этих значений. Есть ли способ извлечь это значение из решателя?

Пример системы приведен ниже: 

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"
...