Prematurely stopping integration in R's deSolve when certain condition is met without using roo-CodePudding

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

CodePudding user response：

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