r - ODE 时代 Matlab 与 R

Question

如果在 matlab 中使用可变时间步长求解器，例如 ODE45 - 我将为输出定义一个时间跨度，即 times = [0 50]，matlab 将返回不同时间步长的结果介于 0 和 50 之间。

但是在 R 中，我似乎必须定义我希望 ODE 返回结果的时间点，即如果我给 times = 0:50，它会在 0,1,2, ... 50。否则，我必须提供一个序列，例如 times = seq(0,50,0.1)。

我有一个功能，它在开始时变化很快，然后逐渐变化。在 MATLAB 中，输出结果反射(reflect)了这一点，结果中返回了 82 个时间步长，其中 49 个在时间步长 0 和 1 之间。

我想知道是否有一种方法可以让 R 以与 MATLAB 相同的方式返回结果，因此无需我预先指定我希望返回结果的时间点。

Answer 1

R-Package deSolve 以这种方式描述了使用的参数times:

time sequence for which output is wanted;

Dennis like to the document还有一个重要的句子:

We note that, for some implementations, the vector times at which the output is wanted defines the mesh at which the method performs its steps, so the accuracy of the solution strongly depends on the input vector times.

下面是一个简单的例子，洛伦兹方程，在关于包 deSolve 的书中提到:

library(deSolve)

parameters <- c(
  a = -8/3,
  b = -10,
  c =  28)

state <- c(
  X = 1,
  Y = 1,
  Z = 1)

# ---- define function in R
Lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    # rate of change
    dX <- a*X + Y*Z
    dY <- b * (Y-Z)
    dZ <- -X*Y + c*Y - Z

    # return the rate of change
    list(c(dX, dY, dZ))
  })   # end with(as.list ...
}

times_1 <- seq(0, 100, by = 1)
out_1 <- lsoda(y = state, times = times_1, func = Lorenz, parms = parameters)

times_2 <- seq(0, 100, by = 0.01)
out_2 <- lsoda(y = state, times = times_2, func = Lorenz, parms = parameters)

tail(out_1)

       time        X          Y          Z
 [96,]   95 30.54833  11.802554  12.445819
 [97,]   96 21.26423   4.341405   4.785116
 [98,]   97 33.05220  13.021730  12.934761
 [99,]   98 21.06394   2.290509   1.717839
[100,]   99 10.34613   1.242556   2.238154
[101,]  100 32.87323 -13.054632 -13.194377

tail(out_2)

           time        X          Y         Z
 [9996,]  99.95 17.16735  -7.509679 -12.08159
 [9997,]  99.96 17.66567  -7.978368 -12.77713
 [9998,]  99.97 18.26620  -8.468668 -13.47134
 [9999,]  99.98 18.97496  -8.977816 -14.15177
[10000,]  99.99 19.79639  -9.501998 -14.80335
[10001,] 100.00 20.73260 -10.036203 -15.40840

您可以在 t = 100 时看到结果的差异。这来自不同定义的时间。

问候，
J_F