我知道 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 后的积分执行了两次。如果满足条件后停止第一次积分,然后直接进行第二次积分,效率会更高。这就是为什么我问是否有任何方法可以在满足特定条件时停止集成,并返回到该点的结果
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
发布评论
评论(1)
首先,deSolve的几个求解器支持根查找,可以在包vignettes或此处。
在其他情况下,始终可以将 ode 嵌入到自己的函数中。这里有一种支持任何求解器的可能方法,它是在OP提供代码示例之前编写的。它使用相当通用的方法,因此应该可以使其适应特定问题。这个想法是在循环内以“走走停停”的模式运行求解器,其中积分的结果用作下一个时间步的初始值。
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.