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

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

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

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

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

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

注意：我的真实方程更复杂，实际上计算起来相当复杂——它们与化学无关，除了值都是正数，一旦接近零，没有赔率，状态系统停止更改。

我有一个手卷的 RK4 实现，可以做我想要的，但我写的代码是我需要测试的代码，我知道 deSolve 包已经过测试，所以我只是在寻找一种方法以固定步长获取输出直到“根”函数返回真

更新

以下是解决我的问题的示例，其中我知道我想要答案的时间：

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,rootfun=function(t,p){y.B<1e-8}
)

其中 __name_of_function__ 希望是 deSolve 包或相关帮助包中的某个函数

解决方法

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,{
    ifelse(B<1e-8,B)
  })
}

# First find where the root is
out <- ode(y0,c(0,1e8),kernel,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,root = rootfun)

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