在不使用根的情况下满足特定条件时过早地停止 R 的 deSolve 中的集成

Question

我知道 R 的 deSolve 中的提前终止可以通过使用根函数而不提供事件函数来实现，这将导致在找到根时终止积分。但是，通过使用此过程，我们仅限于应用具有求根功能的求解器。

事实上，我正在处理一个根的确切位置无关紧要的问题。我需要引入状态变量的突然变化，但发生这种情况的确切时刻并不重要。所以我可以在满足条件时停止积分，重新计算引入突然变化的新起始状态向量，然后重新开始积分。这仍然可以让我灵活地使用 deSolve 包提供的许多解算器中的任何一个。

有推荐的方法吗？

编辑

让我们考虑以下简化示例。所表示的系统是在一维中以 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)

但是，出于上述原因，我正在尝试不使用 root 来执行此操作。可以通过在输出中包含一些变量（必须是数字，否则 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 的几个求解器支持求根，table 可以在包 vignettes 或 here.

中找到

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

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

在不使用根的情况下满足特定条件时过早地停止 R 的 deSolve 中的集成

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

r

ode

differential-equations

desolve