r - deSolve 中 Runge-Kutta 方法 ode45 的自适应时间步长

Question

我想使用 deSolve R 包中的显式 Runge-Kutta 方法 ode45(别名 rk45dp7)来求解具有可变步长的 ODE 问题。

根据 deSolve 文档，可以为 rk 使用自适应或可变时间步长 使用 ode45 方法而不是等距时间步的求解器函数，但我不知道如何执行此操作。

rk 函数是这样调用的:

rk(y, times, func, parms, rtol = 1e-6, atol = 1e-6, verbose = FALSE, tcrit = NULL,
hmin = 0, hmax = NULL, hini = hmax, ynames = TRUE, method = rkMethod("rk45dp7", ... ), 
maxsteps = 5000, dllname = NULL, initfunc = dllname, initpar = parms, rpar = NULL, 
ipar = NULL, nout = 0, outnames = NULL, forcings = NULL, initforc = NULL, fcontrol = 
NULL, events = NULL, ...)

时间是需要明确估计 y 的时间向量。

对于距离为 0.01 的等距时间步，我可以将时间写为

times <- seq(0, 100, 0.01)

假设我想求解 0 到 100 区间的方程，我该如何定义时间而不给出步长？

如有任何帮助，我们将不胜感激。

Answer 1

这里有两个问题。首先，如果您想指定具有多个增量的时间向量，请使用此(例如):

times <- c(seq(0,0.9,0.1),seq(1,10,1))
times
#  [1]  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0 10.0

这里，我们有 [0,1] 乘以 0.1，[1,10] 乘以 1。

但实际上您不必这样做:参数 times=告诉rk(...)报告结果的时间。自适应算法将内部调整时间增量，以在参数指定的时间产生准确的结果。所以对于自适应算法，例如 method="rk45dp7"你什么都不用做。对于非自适应算法，例如 method="euler" ，算法使用的时间增量确实是times=中指定的增量.您可以在这个集成了范德波尔振荡器的简单示例中看到其效果。

y.prime <- function(t,y.vector,b) {    # Van der Pol oscillator
  x <- y.vector[1]
  y <- y.vector[2]
  x.prime <- y
  y.prime <- b*y*(1-x)^2 - x
  return(list(c(x=x.prime,y=y.prime)))
}
h  <- .001                  # time increment
t  <-  seq(0,10,h)          # times to report results
y0 <- c(0.01,0.01)          # initial conditions
euler   <- rk(y0, t,func=y.prime,parms=1,method="euler")
rk45dp7 <- rk(y0, t,func=y.prime,parms=1, method="rk45dp7")
# plot x vs. y
par(mfrow=c(1,2))
plot(euler[,2],euler[,3], type="l",xlab="X",ylab="Y",main=paste("Euler: h =",format(h,digits=3)))
plot(rk45dp7[,2],rk45dp7[,3], type="l",xlab="X",ylab="Y",main=paste("RK45DP7: h =",format(h,digits=3)))

下面是对 h 的几个值的结果比较.请注意 method="rk45dp7" h < 0.5 的结果稳定.这是因为 rk45dp7正在根据需要在内部调整时间增量。对于 method="euler"结果不匹配 rk45dp7直到 h~0.01 .