如何以固定步骤获得输出？

Question

我有一个微分方程系统，它有点复杂，但我能得出的最接近的类比是催化化学反应，其中催化剂随着时间的推移而分解，例如

A + B -> C + B  (rate k1)
B -> D          (rate k2)

所以我们有dAdt(A,B) = -k1*A*B 和dBbt(a,B) = -k2*B

我可以用rk4 和一组固定的时间点来建模...但我想运行这个直到说 B

我认为这可以通过使用通用 ode 包中的根函数来完成，但这也只需要一个 times 参数，这迫使我提前选择我想要“A 和 B 的浓度”的时间.

注意：我的真实方程更复杂，实际上计算起来也很复杂——它们与化学无关，除了值都是正数，一旦一个接近零且没有任何可能性，则系统停止变化。

我有一个手动 RK4 实现，可以满足我的要求，但我编写的代码是我需要测试的代码，我知道 deSolve 包已经过良好测试，所以我只是在寻找一种方法来获得以固定步长输出直到“根”函数返回 true

更新

这是一个解决我的问题的例子，我知道我想要答案的时间：

kernel <- function(t, y, params, input) {
  with(as.list(c(params,y)), {
    dA <- -k1*A*B
    dB <- -k2*B
    list(c(dA,dB))
  })
}
params <- c(k1=0.03, k2=0.005)
y0 <- c(A=1.3, B=0.5)

# here I need to pick the times!
times <- seq(0,500, by=1)
out <- rk4(y0, times, kernel, params)

我想要的是“将代码的最后两行更改为此”

out <- ___name_of_function___(
  y0, 
  initial_time=0, 
  delta_time=1, 
  kernel, 
  params, 
  rootfun=function(t,y,p){y.B<1e-8}
)

__name_of_function__ 希望是 deSolve 包或相关帮助程序包中的某些功能

Answer 1

https://www.rdocumentation.org/packages/deSolve/versions/1.28/topics/lsodar 中的示例 2 表明，一般方法暗示，对于求解器的两次调用，我们可以让第二次调用评估所有需要的点，但代价是无法控制第一次调用中的评估次数：

kernel <- function(t, y, params, input) {
  with(as.list(c(params,y)), {
    dA <- -k1*A*B
    dB <- -k2*B
    list(c(dA,dB))
  })
}
params <- c(k1=0.03, k2=0.005)
y0 <- c(A=1.3, B=0.5)
rootfun <- function(t, y, params, input) {
  with(as.list(c(params,y)), {
    ifelse(B<1e-8,0,B)
  })
}

# First find where the root is
out <- ode(y0, c(0,1e8), kernel, params, root=rootfun)
# Now the second rwo of `out` will be the location of the root
# So just create the times up to that point
times <- seq(0,round(out[2,'time']), by=1)
out <- ode(y0, times, kernel, params, root = rootfun)

注意：求根期间函数的评估次数受 maxsteps 限制，默认为 500