问题描述
我有一个微分方程组,有点复杂,但我能得出的最接近的类比是催化化学反应,其中催化剂随时间分解,例如
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