Jak rozwiązać ODE z progiem wewnętrznym?

Mam następujący funkcja zawierająca pewne ody:Jak rozwiązać ODE z progiem wewnętrznym?

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) 
    }) 
}

Oto niektóre wyniki:

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 
...

Zjednoczone Y są różne, może ktoś mi wyjaśnić, dlaczego proszę?

Uważam, że nie ustawiłem poprawnie progu. Czy jest na to sposób?

Dzięki!

Źródło

2012-10-31 Claudia

Próbowałem odszyfrować plik pomocy dla 'ode' ale bez powodzenia. Wszystko, co mogę zasugerować, to wypróbować bardzo prostą funkcję i zobaczyć, co reprezentują różne zwrócone kolumny. Podejrzewam, że dwie kolumny "Y" patrzą na różne rzeczy. –

Poprawiłem kod, aby podkreślić, że kolumny 3 i 5 powinny wyglądać tak samo. Jednak kiedy próg staje się aktywny (z linii 11), zwracają różne wyniki ?! – Claudia

Ponownie, nie wiem, co to jest podstawowy algorytm "sodowy *", ale domyślam się, że '10590.600567' jest wynikiem poprzedniego cyklu (gdy' dY' wciąż było niezerowe), a ta wartość pozostaje w czymkolwiek kolumna "Y input" reprezentuje, podczas gdy "Y output" jest prawidłowo zamrożone na 35000. –

Pomyśl najprostsza metoda rozwiązywania równań różniczkowych, czyli metoda Eulera:

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

h to mały krok czasowy, myfunction jest f'(t) pochodną f(t) i tylko ocenia pochodnej, nie ma dostępu do rzeczywiste state ani Y. Oba są ustawione wewnętrznie w ode przy użyciu metody, która w zasadzie jest podobna do metody Eulera: biorąc pod uwagę wartości liczbowe f(t),f'(t),h, po prostu aktualizuje się stan f(t+h).

Tak więc próg dostosowuje dY, ale nie może uzyskać dostępu do state["Y"]. Proces prostu manipuluje zmienną lokalną, która jest oceniana jako 35000 w dX <- a*X^6 + Y*Z i dZ <- -X*Y + c*Y - Z ale rzeczywistej state["Y"] są nadpisywane po myfuction powrócił wewnątrz funkcji ode.

Obawiam się, że nie mogę wymyślić prostego sposobu na ominięcie tego projektu. Po prostu użyłbym out[5].

Źródło

2013-07-08 19:49:52 user2561915

Jak rozwiązać ODE z progiem wewnętrznym?

Odpowiedz

Powiązane problemy