在R's deSolve中不使用根的情况下提前停止积分

Question

我知道R的deSolve中的过早终止可以通过使用根函数而不提供事件函数来实现，这将导致在找到根时终止集成。然而，通过使用这一过程，我们仅限于应用具有寻根能力的求解器。

事实上，我所处理的是一个根本位置并不重要的问题。我需要在状态变量中引入一个突然的变化，但是发生这种情况的确切时刻并不重要。因此，我可以在满足条件时停止积分，通过引入突变重新计算新的启动状态向量，然后再开始积分。这仍将使我能够灵活地使用通过deSolve包可用的许多解决程序中的任何一个。

有推荐的方法吗？

编辑

让我们考虑下面的简化示例。所表示的系统是一个1维运动的物体，其速度常数为1，物体从x=0位置开始，并沿维数的正方向运动。我们的目的是改变坐标的原点，例如当物体与原点的距离达到10或更高时，参照x=10的位置，这可以简化为从这个位置减去10。

通过使用根，可以实现以下目标：

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)

然而，出于上述原因，我试图在不使用根的情况下做到这一点。可以通过在输出中包含一些变量(必须是数值变量，因为否则deSolve将不接受它作为每个评估时间的输出)来实现这一点，从而跟踪条件是否已经满足，然后检查集成的结果以确定条件何时满足，应用所需的更改，从该点重新执行集成，然后组合输出。使用与前面相同的示例：

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)

但是，这会导致无用的计算，因为在time=10之后执行两次集成。如果在满足条件后停止第一次集成，然后直接进行到第二集成部分，则会更有效。这就是为什么我要问，在满足某一条件时，是否有任何方法来停止集成，并将结果返回到这一点。

Answer 1

首先，几个deSolve的求解者支持根查找，一个表可以在包中找到小片段或这里。

在其他情况下，总是可以将ode嵌入到自己的函数中。这里，一种支持任何求解器的可能方法，是在OP之前编写的，它提供了一个代码示例。它使用了一种相当通用的方法，因此应该可以根据具体的问题对其进行调整。其思想是在循环中以“停止和走”模式运行求解器，其中集成的结果被用作下一个时间步骤的初始值。

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