Wie löst man ODEs mit einem internen Schwellenwert?

Ich habe die folgende Funktion einige Oden enthält:Wie löst man ODEs mit einem internen Schwellenwert?

myfunction <- function(t, state, parameters) { 
    with(as.list(c(state, parameters)),{ 
     if (X>20) {   # this is an internal threshold! 
      Y <- 35000 
      dY <- 0 
     }else{ 
      dY <- b * (Y-Z) 
     } 
     dX <- a*X^6 + Y*Z 
     dZ <- -X*Y + c*Y - Z 
     # return the rate of change 
     list(c(dX, dY, dZ),Y,dY) 
    }) 
}

Hier einige Ergebnisse:

library(deSolve) 
parameters <- c(a = -8/3, b = -10, c = 28) 
state <- c(X = 1, Y = 1, Z = 1) 
times <- seq(0, 10, by = 0.1) 
out <- ode(y = state, times = times, func = myfunction, parms = parameters) 
out 
    time   X   Y   Z   Y   dY 
1 0.0 1.000000  1.000000  1.000000  1.000000  0.00000 
2 0.1 1.104670  2.132728  4.470145  2.132728 23.37417 
3 0.2 1.783117  6.598806 14.086158  6.598806 74.87351 
4 0.3 2.620428 20.325966 42.957134 20.325966 226.31169 
5 0.4 3.775597 60.969424 126.920014 60.969424 659.50590 
6 0.5 5.358823 176.094907 358.726482 176.094907 1826.31575 
7 0.6 7.460841 482.506706 953.270570 482.506706 4707.63864 
8 0.7 10.122371 1230.831764 2330.599161 1230.831764 10997.67398 
9 0.8 13.279052 2859.284114 5113.458479 2859.284114 22541.74365 
10 0.9 16.711405 5912.675147 9823.406760 5912.675147 39107.31613 
11 1.0 24.452867 10590.600567 16288.435139 35000.000000  0.00000 
12 1.1 25.988924 10590.600567 23476.343542 35000.000000  0.00000 
13 1.2 26.572411 10590.600567 26821.703961 35000.000000  0.00000 
14 1.3 26.844240 10590.600567 28510.668725 35000.000000  0.00000 
15 1.4 26.980647 10590.600567 29391.032472 35000.000000  0.00000 
...

Staaten Y verschieden sind, kann mir jemand erklären, warum bitte?

Ich glaube, ich habe meine Schwelle nicht korrekt eingestellt. Gibt es einen Weg dazu?

Danke!

Quelle

2012-10-31 Claudia

Ich habe versucht, die Hilfedatei für 'Ode' zu entschlüsseln, aber ohne Erfolg. Alles, was ich vorschlagen kann, ist eine sehr einfache Funktion zu versuchen und zu sehen, was die verschiedenen zurückgegebenen Spalten darstellen. Ich vermute, dass die zwei "Y" -Säulen verschiedene Dinge betrachten. –

Ich habe den Code geändert, um hervorzuheben, dass die Spalten 3 und 5 beide dasselbe betrachten sollten. Wenn der Schwellenwert jedoch aktiv wird (ab Zeile 11), liefern sie unterschiedliche Ergebnisse ?! – Claudia

Noch einmal, ich weiß nicht, was der zugrundeliegende 'soda *' Algorithmus tut, aber meine Vermutung ist nun, dass '10590.600567' die Ausgabe des vorherigen Zyklus ist (wenn' dY' immernoch nicht Null war), und dieser Wert bleibt in was auch immer Die Spalte "Y-Eingabe" stellt dar, während die "Y-Ausgabe" ordnungsgemäß bei 35000 eingefroren ist. –

Denken die einfachste Methode ODEs, dh Euler-Verfahren zu lösen:

state = state+myfunction(t,state,parameters)*h 
f(t+h)=f(t) + f'(t) *h

h ist ein kleiner Zeitschritt myfunction das f'(t) Derivat von f(t) ist und wertet nur die Ableitung, es hat keinen Zugriff auf die tatsächliche state noch Y. Beide werden intern in ode unter Verwendung einer Methode eingestellt, die im Prinzip ähnlich zu Euler ist: bei den numerischen Werten von f(t),f'(t),h aktualisiert sie einfach den Zustand f(t+h).

Der Schwellenwert passt also dY an, kann aber nicht auf state["Y"] zugreifen. Der Prozess manipuliert nur eine lokale Variable, die als 35000 in dX <- a*X^6 + Y*Z und dZ <- -X*Y + c*Y - Z ausgewertet wird, aber die tatsächliche state["Y"] wird überschrieben, nachdem die myfuction innerhalb der ode Funktion zurückgekehrt ist.

Ich fürchte, dass ich nicht an eine einfache Möglichkeit denke, dieses Design zu umgehen. Ich würde einfach out[5] verwenden.

Quelle

2013-07-08 19:49:52 user2561915

Wie löst man ODEs mit einem internen Schwellenwert?

Antwort

Verwandte Themen