deSolve in R дает неправильные промежуточные значения

Question

У меня возникают проблемы с решателем ode в R от deSolve.

Я обнаружил, что решатель дает ненулевые промежуточные значения для переменных, которые всегда должны быть равны нулю. Это не похоже на проблему, пока я не начну отлаживать тесты, например if (R> 0) {browser()}, которые запускаются.

Мой код ниже. Заранее спасибо!

Эллен

library(deSolve) 

simpleSIR <- function(t,states,par){ 
    with(as.list(c(states,par)),{ 

    S=states[1] 
    I=states[2] 
    R=states[3] 

    newinfections = beta*I*S 

    dS <- -newinfections 
    dI <- newinfections - gamma*I 
    dR <- gamma*I 

    if(R>0) 
    { 
     print(paste("why is this not zero?",R)) 
    } 

    return(list(c(dS,dI,dR))) 
})} 


par=list(beta=0.3,gamma=0.0) 

init=c(0.99,0.01,0) 

times <- seq(0,500,by = 1) 


out <- as.data.frame(ode(y = init, times = times, func = simpleSIR, parms = par,maxsteps=2000)) 


plot(out[,2],type="l") 
lines(out[,3],type="l",col="red") 
lines(out[,4],type="l",col="blue") 

Answer 1

Я думаю, ваша проблема заключается в том, что метод по умолчанию для численного интегрирования в ode() является lsoda. Этот метод может переключаться между решателями для жестких и нежестких проблем. В случае, если вы переключаетесь на жесткие решатели, якобиан вычисляется численно, что может привести к появлению числовых ошибок. Вы можете видеть, что это может быть причиной, в следующем коде:

out <- deSolve::ode(y = init, times = times, func = simpleSIR, parms = par,maxsteps=2000) 
deSolve::diagnostics.deSolve(out) 

«[...] 14 Количество оценок якобиевыми и LU разбиений до сих пор: 23 [...]»

, который соответствует количеству печатных сообщений (23), которые производит ваш исходный код. Вы избавитесь от проблемы с помощью нежестких решателей как RK4:

out.rk4 <- deSolve::ode(y = init, times = times, func = simpleSIR,method = "rk4", parms = par,maxsteps=2000) 

, если вы настаиваете на использование lsoda, вы можете попробовать поставлять lsoda с якобианом Вы вычислили аналитически.

Answer 2

Прежде всего, я не знаю, если вы понимаете различия между государственными-переменными и параметрами.

R (в вашем случае) является переменной состояния. Это означает, что в начале моделирования (при t = 0) R = 0. Но с развитием времени это может измениться. См. Переменную состояния S!

Ваши параметры бета и гамма, они не будут меняться во время вычислений.

Тогда вы можете упростить свой код. Используйте эти строки вместо вашего определения par и init:

par <- c(beta = 0.3, gamma = 0.0) 

init <- c(S = 0.99, I = 0.01, R = 0) 

И теперь вы можете удалить следующие строки, потому что с with(as.list(c(states, par)) имена параметров и государственных параметров доступны с их именами S, I и R .

S = states[1] 
I = states[2] 
R = states[3] 

И вы можете упростить сюжет, тоже:

matplot(out[,1], out[,2:4], type = "l", col = c("black", "red", "blue")) 

С уважением,

J_F