Làm thế nào để giải quyết ODE với ngưỡng nội bộ?

Question

Tôi đã hàm sau có chứa một số Odes:

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

Dưới đây là một số kết quả:

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

Hoa Y là khác nhau, có thể ai giải thích cho tôi tại sao không?

Tôi tin rằng tôi chưa đặt đúng ngưỡng của mình. Có cách nào không?

Cảm ơn!

Answer 1

suy nghĩ phương pháp đơn giản nhất để giải quyết ODEs, tức là phương pháp Euler:

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

h là một bước thời gian nhỏ, myfunction là f'(t) phái sinh của f(t) và chỉ đánh giá đạo hàm, nó không có quyền truy cập vào các thực tế state cũng không phải Y. Cả hai được thiết lập nội bộ trong ode sử dụng một phương pháp mà về nguyên tắc tương tự như Euler: cho các giá trị số của f(t),f'(t),h nó chỉ cập nhật trạng thái f(t+h).

Do đó ngưỡng điều chỉnh dY nhưng không thể truy cập state["Y"]. Quá trình này chỉ thao tác biến cục bộ được đánh giá là 35000 trong dX <- a*X^6 + Y*Z và dZ <- -X*Y + c*Y - Z nhưng thực tế state["Y"] bị ghi đè sau khi myfuction đã trở lại bên trong hàm ode.

Tôi sợ rằng tôi không thể nghĩ ra một cách đơn giản để bỏ qua thiết kế này. Tôi chỉ sử dụng out[5].