3

我想使用deSolve R 包中的显式 Runge-Kutta 方法ode45（别名 rk45dp7）来解决可变步长的 ODE 问题。

根据 deSolve 文档，可以使用 ode45 方法对rk求解器函数使用自适应或可变时间步长，而不是等距时间步长，但我不知道如何做到这一点。

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, ...)

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

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

times <- seq(0, 100, 0.01)

假设我想求解从 0 到 100 的区间的方程，如何在不给出步长的情况下定义时间？

任何帮助将不胜感激。

4

1 回答 1

5

这里有两个问题。首先，如果要指定具有多个增量的时间向量，请使用此（例如）：

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] x 0.1 和 [1,10] x 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。

于 2013-12-15T16:12:47.910 回答