Jacobian für verschiedene Variablen

Question

eine Differentialgleichung könnte als

sys <- function(t, y, p, u) { 
    dy <- numeric(2) 
    u <- u(t) 

    dy[1] = p$a*(1 + p$b*(y[2] - 1)/(p$c + y[2] - 1) - u*y[1]) 
    dy[2] = u*y[1] - y[2] 
    list(dy) 
} 

Lassen Sie uns definiert werden Weiteren davon aus, dass die stationären Zustände (Gleichgewicht) sind bekannt. Nun, gibt es eine Möglichkeit, die Jacobi-Matrix der rechten Seite in Bezug auf x zu berechnen?

Ich weiß, dass ich

f <- function(y){ 
    c(
    p$a*(1 + p$b*(y[2] - 1)/(p$c + y[2] - 1) - u*y[1]), 
    u*y[1] - y[2] 
    ) 
} 

und berechnen Sie die Jacobi mit

Jx <- jacobian(f, c(1,1)) 

definieren könnte, wobei jacobian aus dem pracma Paket kommt. Aber gibt es keinen leichteren Weg ohne diesen Zwischenschritt? Es könnte auch helfen, wenn man f innerhalb sys, z.B.

sys <- function(t, y, p, u) { 
    dy <- numeric(2) 
    u <- u(t) 

    dy[1] = f(y)[1] 
    dy[2] = f(y)[2] 
    list(dy) 
} 

Und zuletzt, könnte es auch eine Möglichkeit geben, die Jacobi w.r.t zu berechnen. zu u?

Vielen Dank!

Answer 1

Es gibt eine Lösung mit dem R Paket rootSolve.

Dafür muss Ihre Funktionsdefinition ein wenig anders sein, aber (meiner Meinung nach) bequemer. Ich weiß nicht, Ihre genauen Parameter p oder Ihre Funktion u, also machte ich ein minimales Beispiel:

library(rootSolve) 
sys <- function(t, y, parms) { 
    with(as.list(c(y,parms)),{ 

    dy = a*(1 + b*(z - 1)/(c + z - 1) - 1*y) 
    dz = 1*y - z 
    return(list(c(dy, dz))) 
    }) 
} 

parms <- list(a = 1, b = 1, c= 2) 

rootSolve::jacobian.full(y = c(y = 1, z = 1), func = sys, parms = parms) 

In der Funktion jacobian.fully() Sie können Ihre steady state Ergebnisse verwenden, nahm ich nur zufällige Ergebnisse. Die Definition von sys verwendet die Standarddefinition für ODEs im Paket deSolve, ein ausgezeichnetes Paket zur Lösung von ODEs.
Das Ergebnis ist eine normale Jacobi-Matrix.

Also mit dieser Definition können Sie lösen Algorithmen für Ihre Gleichungen, z.

library(deSolve) 
ode <- deSolve::ode(y = c(y = 1, z = 0), 
      times = seq(1,100), 
      func = sys, 
      parms = parms) 
plot(ode) 

Ich hoffe, das hilft Ihnen ein wenig!
Grüße,
J_F