ODE 时代 Matlab 与 R

Question

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

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

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

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

Answer 1

在 deSolve 上阅读 this document 后，它指出：

We solve the IVP for 100 days, producing output every 0.01 days; this small output step is necessary to obtain smooth figures. In general this does not affect the time step of integration; this is usually determined by the solver

因此，您似乎确实应该使用 times = seq(0,50,0.1) 作为输入。如果您只想在图表中显示 'interesting' 个点，我想您可以编写一个小函数来 post-处理实际求解器输出的输出。

Answer 2

R-Package deSolve这样描述使用的参数times：

time sequence for which output is wanted;

Dennis 喜欢的文档还有一个重要的句子：

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 时看到结果的差异。这来自不同的定义 times。

此致，
J_F

ODE 时代 Matlab 与 R

ODE Times Matlab vs R

matlab

r

ode45