3

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

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

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

我想知道是否有办法让 R 以与 MATLAB 相同的方式返回结果，所以没有我预先指定我希望返回结果的时间点。

4

2 回答 2

2

R -Package以这种方式deSolve描述使用的参数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 时看到结果的差异。这来自不同的定义times。

问候，
J_F

于 2016-06-20T09:57:56.570 回答
1

在deSolve 上阅读此文档后，它指出：

我们解决IVP 100天，每0.01天产出；这个小的输出步骤对于获得平滑的数字是必要的。一般来说，这不会影响积分的时间步长；这通常由求解器确定

因此，您似乎确实应该将其times = seq(0,50,0.1)用作输入。如果您只想在图表中显示“有趣”的点，我想您可以编写一个小函数来对实际求解器输出的输出进行后处理。

于 2016-06-20T09:09:43.833 回答