Ist es möglich, rk4 und rootfun in ode (Paket deSolve) zu verwenden

Question

Ich versuche, ein Modell von Beute-Beute-Raubtier-System mit Differentialgleichungen basierend auf dem LV-Modell zu modellieren. Aus Gründen der Präzision muss ich die runge-kutta4-Methode verwenden. Aber angesichts der Gleichungen werden einige der Populationen schnell negativ. So

Ich versuchte, die Ereignisse/Wurzelsystem von ODE zu verwenden, aber es scheint, dass rk4 und rootfun sind nicht kompatible Geräte ...

eventFunc <- function(t, y, p){ 
    if (y["N1"] < 0) { y["N1"] = 0 } 
    if (y["N2"] < 0) { y["N2"] = 0 } 
    if (y["P"] < 0) { y["P"] = 0 } 
    return(y) 
} 

rootFunction <- function(t, y, p){ 
    if (y["P"] < 0) {y["P"] = 0} 
    if (y["N1"] < 0) {y["N1"] = 0} 
    if (y["N2"] < 0) {y["N2"] = 0} 
    return(y) 
} 

out <- ode(func=Model_T2.2, 
      method="rk4", 
      y=state, 
      parms=parameters, 
      times=times, 
      events = list(func = eventFunc, 
         root = TRUE), 
      rootfun = rootFunction 
      ) 

Dieser Code gibt mir die followin Fehler:

Fehler in checkevents (events, times, Yname, dllname): entweder 'events $ time' sollte angegeben werden und die Zeiten der Events enthalten, wenn 'events $ func' angegeben ist und keine Root-Funktion oder Ihr Solver root nicht unterstützt Funktionen

Gibt es eine Lösung, um rk4 zu verwenden und zu verbieten, dass die Funktionen unter 0 gehen?

Vielen Dank im Voraus.

Für diejenigen, die fragen könnte, hier ist das, was funktioniert:

if(!require(ggplot2)) { 
    install.packages("ggplot2"); require(ggplot2)} 
if(!require(deSolve)) { 
    install.packages("deSolve"); require(deSolve)} 
Model_T2.2 <- function(t, state, par){ 
    with(as.list(c(state, par)), { 

    response1 <- (a1 * N1)/(1+(a1*h1*N1)+(a2*h2*N2)) 
    response2 <- (a2 * N2)/(1+(a1*h1*N1)+(a2*h2*N2)) 

    dN1 = r1*N1 * (1 - ((N1 + A12 * N2)/K1)) - response1 * P 
    dN2 = r2*N2 * (1 - ((N1 + A21 * N2)/K2)) - response2 * P 
    dP = ((E1 * response1) + (E2 * response2)) * P - Mp 

    return(list(c(dN1, dN2, dP))) 
    }) 
} 

parameters<-c(
    r1=1.42, r2=0.9, 
    A12=0.6, A21=0.5, 
    K1=50, K2=50, 
    a1=0.77, a2=0.77, 
    b1 = 1, b2=1, 
    h1=1.04, h2=1.04, 
    o1=0, o2=0, 
    Mp=0.22, 
    E1=0.36, E2=0.36 
) 

## inital states 
state<-c(
    P=10, 
    N1=30, 
    N2=30 
) 

times <- seq(0, 30, by=0.5) 

out <- ode(func=Model_T2.2, 
      method="rk4", 
      y=state, 
      parms=parameters, 
      times=times, 
      events = list(func = eventFunc, 
         root = TRUE), 
      rootfun = rootFunction 
      ) 

md <- melt(as.data.frame(out), id.vars=1, measure.vars = c("N1", "N2", "P")) 
pl <- ggplot(md, aes(x=time, y=value, colour=variable)) 
pl <- pl + geom_line() + geom_point() + scale_color_discrete(name="Population") 
pl 

Und das Ergebnis in einem Diagramm: Evolution of prey1, prey2 and predator populations Wie Sie sehen können, die Bevölkerung von Räubern negativ werden, die in der realen offensichtlich unmöglich ist Welt.

Edit: fehlende Variablen, tut mir leid.

Answer 1

Dies ist ein Problem, das Sie mit allen expliziten Lösern wie rk4 haben werden. Die Reduzierung des Zeitschritts hilft bis zu einem gewissen Punkt. Verwenden Sie besser einen Solver mit einer impliziten Methode, lsoda scheint universell in der einen oder anderen Form verfügbar.

Eine weitere Möglichkeit, positive Werte explizit zu erzwingen, besteht darin, sie als Exponentialwerte zu parametrisieren. Set N1=exp(U1), N2=exp(U2) übersetzt dann die ODE Funktionscode (als dN = exp(U)*dU = N*dU)

N1 <- exp(U1) 
N2 <- exp(U2) 
response1 <- (a1)/(1+(a1*h1*N1)+(a2*h2*N2)) 
response2 <- (a2)/(1+(a1*h1*N1)+(a2*h2*N2)) 

dU1 = r1 * (1 - ((N1 + A12 * N2)/K1)) - response1 * P 
dU2 = r2 * (1 - ((N1 + A21 * N2)/K2)) - response2 * P 
dP = ((E1 * response1*N1) + (E2 * response2*N2)) * P - Mp 

Für die Ausgabe, die Sie dann natürlich haben U1, U2N1, N2 aus den Lösungen zu rekonstruieren.

Answer 2

Dank J_F kann ich jetzt mein L-V Modell laufen lassen.

Die radau (nicht randau wie Sie erwähnt) Funktion tatsächlich akzeptieren Wurzelfunktion und Ereignisse ans implementiert implizit die Runge-Kutta-Methode.

Nochmals vielen Dank, hoffe, dass dies jemand in der Zukunft helfen wird.