Prematurely stopping integration in R's deSolve when certain condition is met without using roots

Question

I know that premature termination in R's deSolve can be achieved by using root functions and not providing an event function, which will result in termination of integration when a root is found. However, by using this procedure, we are limited to applying a solver with root-finding abilities.

I am in fact dealing with a problem for which the exact position of the root does not matter. I need to introduce a sudden change in the state variables, but the exact moment when this happens is not important. So I can just stop the integration when the condition is met, recalculate a new starting state vector with the sudden change introduced, and start integration again. This would still give me the flexibility of being able to use any of the many solvers available through the deSolve package.

Is there a recommended way to do this?

Edit

Let's consider the following simplified example. The represented system is an object moving in 1 dimension with constant velocity of 1. The object starts at position x=0, and moves in the positive direction of the dimension. We aim to perform a change of the origin of coordinates such than when the object reaches a distance of 10 or higher from the origin, the position is referenced with respect to the point at which x=10. This can be simplified as subtracting 10 from the position at this point.

Using roots, this can be achieved as follows:

library(deSolve)
odeModel <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
        x <- state
        dx <- 1
        list(
            dx
        )
    })
}

initialState <- 0
times <- 0:15

rootFunc <- function(t, state, parameters) {
    return(state-10)
}

eventFunc <- function(t, state, parameters) {
    return(state - 10)
}

integrationTest <- ode(y=initialState, times=times, func=odeModel, parms=NULL, 
                       events=list(func=eventFunc, root=TRUE), rootfun=rootFunc)

However, I am trying to do this without using roots for the reasons explained above. It is possible to do it by including in the output some variable (must be numeric, since otherwise deSolve will not accept it as an output of each evaluation time) that keeps track of whether the condition has been fulfilled, then examining the results of integration to identify when the condition was fulfilled, applying the desired change, re-doing integration from that point and then combining the outputs. Using the same example as before:

library(deSolve)
distanceHigherThan10 <- function(x) {
    result <- if(x >= 10) 1 else 0
    return(result)
}

odeModel <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
        x <- state
        dx <- 1
        higherThanThreshold <- distanceHigherThan10(x)
        list(
            dx, higherThanThreshold
        )
    })
}

initialState <- 0
times <- 0:15

integration <- ode(y=initialState, times=times, func=odeModel, parms=NULL)

changePoint <- which(diff(integration[,3]) != 0)[1] + 1
newIntegrationTimes <- times[changePoint:length(times)] - times[changePoint]
newStartingPosition <- integration[changePoint, 2] - 10
newIntegration <- ode(y=newStartingPosition, times=newIntegrationTimes, func=odeModel, parms=NULL)

newIntegration[, 1] <- newIntegration[, 1] + times[changePoint]
newIntegration[, 3] <- 1

combinedResults <- rbind(integration[1:(changePoint-1), ], newIntegration)

However, this leads to useless calculations, since integration after time=10 is performed twice. It would be more efficient to just stop the first integration after the condition has been met, and then proceeding directly to second integration part. This is why I am asking if there is any way to stop the integration when a certain condition is met, returning the results up to that point

It's easier to help you if you include a simple [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input and desired output that can be used to test and verify possible solutions. It's hard to answer questions abstractly and much easier when you actually provide an example we can test with. — MrFlick, Mar 01 '22 at 05:54
The root function should change sign at the intended root. It is doubtful that with the absolute value a root can be found at all, since the solver looks for sign changes over the current integration step. — Lutz Lehmann, Mar 01 '22 at 08:04
@LutzLehmann thanks for pointing this out! Indeed, in my original application, I am not applying absolute values (somehow, the simplified example with absolute value worked nevertheless!). I have edited to correct this — Rafa, Mar 01 '22 at 08:20

tpetzoldt · Accepted Answer · 2022-03-01T07:58:17.000

1

First, several solvers of deSolve support root finding, a table can be found in the package vignettes or here.

In other cases, it is always possible to embed ode in an own function. Here a possible approach that supports any solver, that was written before the OP provided a code example. It uses a rather generic approach, so it should be possible to adapt it to the specific problem. The idea is to run the solver in a "stop-and-go" mode within a loop, where the outcome of the integration is used as initial value of the next time step.

library(deSolve)

## a model with two states
model <- function(t, y, p) {
  list(c(p[1] * y[1] - p[2], p[1] * y[2]))
}

times <- 0:10
p <- c(-0.1, 0.1)
y0 <- c(y1 = 1, y2 = 2)

## standard approach without root detection
out1 <- ode(y = y0, times = times, model, p = p)
out1

## stepper function that stops after root finding
stepping <- function(y0, times, model, p, ...) {
  for (i in 2:length(times)) {
    o <- ode(y = y0, times[c(i-1, i)], func = model, p = p, ...)
    if (i == 2) {
      out <- o[1,]
    }
    ## condition may be adapted
    if (any(o[1, -1] * o[2, -1] < 0)) break 
    y0 <- o[2, -1]
    out <- rbind(out, o[2,])
  }
  out
}

out2 <- stepping(y0, times, model, p)

out1 # all time steps
out2 # stopped at root

edited Mar 01 '22 at 07:58

answered Mar 01 '22 at 07:27

tpetzoldt

5,338
2
12
29

Thanks! I see that the loop is making successive calls to function ode, where at each call only 2 times are passed. Would this affect the performance of solvers that use variable, adaptive step size? – Rafa Mar 01 '22 at 07:35
1

Yes, the repeated re-initialization of `ode` costs performance, where the relative cost depends on size and complexity of the model. Ideed, we use the "stop-and-go"-mode sometimes for complex models, where its cost is negligible. On the other hand, default internal integration steps are also restricted by the steps given in `times`. This design decision was made in the interest of robustness (i.e. to avoid naive mistakes). It can be unlocked by setting `hmax`. – tpetzoldt Mar 01 '22 at 07:56
2

@Rafa : Changing the state in the event function essentially is also a (soft) restart of the integration, as the current step size may be no longer optimal after the state change, the internal state of the step-size controller (if there is some step memory involved) needs to be re-initialized. – Lutz Lehmann Mar 01 '22 at 08:08
I see! My concern was more along the lines of the following. Let's say we have target output times of seq(0,15, by=1). The sudden change in state is introduced at t=10. So for the full blocks before and after that, they can be integrated as a whole. Therefore, integrators with adaptive step size might find the optimal internal step quickly and then keep applying it for the rest of the integration block. Instead, with the proposed "stop-and-go" models, we are only passing 2 target times, and I wonder if this might cause the integrator to never find the optimal internal step size? – Rafa Mar 01 '22 at 08:23
1

I agree with @Lutz Lehmann and want to add that the repeated evaluation of event and root functions have also a cost. It depends on the model and the particular setting what is faster: internal root finding or external stop and go. – tpetzoldt Mar 01 '22 at 08:25
1

Regarding to the comment of @Rafa: The main cost of stop-and-go is not a result of limitation between two time steps. This is the default anyway. The main cost occurs due repeated re-initialization of the solver. My recommendation: try it out and don't worry too much about performance - except if you intend to run the model very often like in MCMC. – tpetzoldt Mar 01 '22 at 08:31
@tpetzoldt Sorry, my bad use of the terminology made me not explain what I meant properly. When I mentioned performance, I meant not computational time/cost, but rather if the solvers will be able to find the optimal internal step size even when provided with just 2 target times. But this is a very interesting point about the performance anyway; I am indeed in a situation where the model will be run many times (calculation of satellite trajectories for long periods). I will compare the costs of "stop-and-go" vs the scheme I described in the example! – Rafa Mar 01 '22 at 08:35
1

The default stepping is always limited between 2 output time steps and then a restart (as far as I remember ;-). This setting was made in of the early days of deSolve to avoid pitfalls and misunderstanding if a solver jumps over the pre-defined steps, while the output values are interpolated, especially for non-autonomous models. It turned out that this works well in practice if the number of output steps is not too large. Experienced users may consider to tweak solver settings for maximum performance. – tpetzoldt Mar 01 '22 at 08:47
@tpetzoldt Thanks a lot for the very insightful comments into the inner workings! These are being extremely helpful. So, if I understood correctly, for example, for the RADAU solver, which has adaptive steps, if we provide as target times c(0,1,2) and c(0,1), the same internal steps between targets t=0 and t=1 are expected to be used in both cases, right? – Rafa Mar 01 '22 at 09:05
Not exactly, e.g. due to rounding errors ... put a `cat(t, "\n")` in the model. – tpetzoldt Mar 01 '22 at 10:27

Prematurely stopping integration in R's deSolve when certain condition is met without using roots

1 Answers1