ODE Times Matlab против R

При использовании решателя с переменным временным шагом, такого как ODE45 в Matlab, я бы определил временной интервал для выходов, то есть times = [0 50], и Matlab будет возвращать результаты с различными временными шагами от 0 до 50.

Однако в R, похоже, мне нужно определить моменты времени, в которые я хочу, чтобы ODE возвращало результаты, т.е. если бы я дал times = 0:50, он вернул бы 51 результат в 0,1,2, ... 50. В противном случае я должен предоставить такую ​​последовательность, как, times = seq(0,50,0.1).

У меня есть функция, которая вначале меняется быстро, а затем гораздо более постепенно. В MATLAB выходные результаты отражают это с помощью 82 временных шагов, возвращаемых в результатах, из которых 49 находятся между временным шагом 0 и 1.

Я хочу знать, есть ли способ заставить R возвращать результаты таким же образом, как и MATLAB, поэтому без предварительного указания моментов времени, в которые я хочу, чтобы результаты возвращались.


r matlab ode45
person Sarah Halliday    schedule 20.06.2016    source источник

Ответы (2)


arrow_upward
2
arrow_downward

Пакет R 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

person J_F    schedule 20.06.2016
comment
Есть ли способ определить лучший временной интервал для аргумента раз? - person Yan; 02.08.2019

arrow_upward
1
arrow_downward

После прочтения этого документ на deSolve гласит:

Мы решаем IVP за 100 дней, производя вывод каждые 0,01 дня; этот небольшой шаг вывода необходим для получения плавных фигур. Обычно это не влияет на временной шаг интегрирования; это обычно определяется решателем

Таким образом, может показаться, что вы действительно должны использовать times = seq(0,50,0.1) в качестве входных данных. Если вы хотите показать только «интересные» точки на графике, я полагаю, вы могли бы написать небольшую функцию для пост-обработки вывода фактического вывода решателя.

person Dennis Jaheruddin    schedule 20.06.2016