deSolve in R gibt falsche Zwischenwerte

Question

Ich habe Probleme mit dem Ode Solver in R von deSolve.

Ich habe gefunden, dass der Löser Zwischenwerte ungleich Null für Variablen gibt, die immer gleich Null sein sollten. Es sieht nicht wie ein Problem aus, bis ich versuche, mit Tests wie if (R> 0) {browser()} zu debuggen, das ausgelöst wird.

Mein Code ist unten. Danke im Voraus!

Ellen

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

Ich denke, das Problem liegt in der Tatsache, dass die Standardmethode für die numerische Integration in ode()lsoda ist. Diese Methode kann zwischen Lösern für steife und nicht steife Probleme wechseln. Falls Sie zu den steifen Solvern wechseln, wird das Jacobi numerisch ausgewertet, was zu den numerischen Fehlern führen kann, die Sie sehen. Sie können sehen, dass dies der Grund in dem folgenden Code sein könnte:

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

"[...] 14 Die Anzahl der Jacobi-Auswertungen und LU Zersetzungen bisher: 23 [...]"

entspricht der Anzahl der Druckmeldungen (23), die Ihr Originalcode erzeugt. Sie loszuwerden, das Problem umgehen, indem ein nicht-steifen Solver wie RK4 mit:

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

wenn Sie sich mit lsoda bestehen, möchten Sie vielleicht Versorgung lsoda mit dem jacobian versuchen Sie analytisch berechnet.

Answer 2

Zunächst einmal, ich weiß nicht, ob Sie die Unterschiede zwischen den Zustandsvariablen und Parameter zu verstehen.

R (in Ihrem Fall) ist eine Statusvariable. Das bedeutet, dass zu Beginn der Simulation (bei t = 0) R = 0. Aber mit der Zeitentwicklung kann sich das ändern. Sehen Sie Ihre Zustandsvariable S!

Ihre Parameter Beta und Gamma ändern sich während der Berechnungen nicht.

Dann können Sie Ihren Code vereinfachen. Verwenden Sie diese Zeilen anstelle von Ihrer Definition von par und init:

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

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

Und jetzt können Sie die folgenden Zeilen entfernen, weil mit with(as.list(c(states, par)) die Namen der Parameter und state-Parameter sind mit ihrem Namen S, I und R .

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

Und Sie können Ihr Grundstück, vereinfachen:

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

Grüße,

J_F