How to solve ODEs 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 * (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)
    })
}

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

States Y are different, can anybody explain me why please?

I believe I haven't set my threshold correctly. Is there a way to that?

Thanks!

Scarce

I tried to decipher the help file for ode but w/o success. All I can suggest is to try a very simple function and see what the various returned columns represent. I suspect the two "Y" columns are looking at different things. – Robet

I have amended the code to highlight that columns 3 and 5 should both look at the same thing. However when the threshold becomes active (from line number 11) they return different results?! – Scarce

Again, I don't know what the underlying soda* algorithm does, but my guess now is that 10590.600567 is the output of the previous cycle (when dY was still nonzero), and that value remains in whatever the "Y input" column represents, while the "Y output" is properly frozen at 35000 . – Robet

?events could be a way to get around this – Doctrine

Answer 1

Think the simplest method to solve ODEs, i.e. Euler 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 f'(t) derivative of f(t) and only evaluates the derivative, it does not have access to the actual state nor Y. Both are set internally in ode using a method which in principle is similar to Euler's: given the numerical values of f(t),f'(t),h it just updates state f(t+h).

So the threshold adjusts dY but cannot access state["Y"]. The process just manipulates a local variable which is evaluated as 35000 in dX <- a*X^6 + Y*Z and dZ <- -X*Y + c*Y - Z but the actual state["Y"] is overwritten after the myfuction has returned inside the ode function.

I am afraid that I cannot think of a simple way to bypass this design. I would just use out[5].

Staggard