1

我知道 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 之后的积分会执行两次。在满足条件后停止第一次积分，然后直接进行第二次积分会更有效。这就是为什么我要问是否有任何方法可以在满足某个条件时停止集成，并将结果返回到该点

4

1 回答 1

1

首先，deSolve的几个求解器支持求根，可以在包 vignettes 或这里找到一个表。

在其他情况下，始终可以嵌入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
于 2022-03-01T07:27:57.737 回答