2

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

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

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

我可以用rk4一组固定的时间点对此进行建模......但我想运行它直到说 B < 1e-8

我认为这可以通过使用通用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或相关的帮助程序包在哪里

4

1 回答 1

0

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

于 2021-04-27T10:01:57.410 回答