How to solve ODE with an internal threshold?

Question

I have the following function containing some 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 * (YZ) } dX <- a*X^6 + Y*Z dZ <- -X*Y + c*Y - Z # return the rate of change list(c(dX, dY, dZ),Y,dY) }) } 

Here are some results:

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

Y states are different, can someone explain to me why?

I believe that I did not set my threshold correctly. Is there any way?

Thanks!

Answer 1

Think of the simplest way to solve ODE, i.e. Euler's method:

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

h is a small time step, myfunction is the derivative of f'(t) from f(t) and only evaluates the derivative, it does not have access to the actual state and Y Both are installed inside ode using a method that is basically similar to Euler's: given the numerical values ​​of f(t),f'(t),h , it just updates the state of f(t+h) .

Thus, the threshold configures dY , but cannot access state["Y"] . The process simply processes the local variable, which evaluates to 35000 in dX <- a*X^6 + Y*Z and dZ <- -X*Y + c*Y - Z , but the actual state["Y"] overwritten after myfuction inside the ode function.

I am afraid that I can’t come up with a simple way around this design. I would just use out[5] .