ODE Times Matlab vs R

Question

If using a variable time step solver such as ODE45 in matlab - I would define a time span for the outputs, i.e. times = [0 50], and matlab would return results at various time steps between 0 and 50.

However in R it appears I have to define the time points at which I want the ODE to return the results, i.e if I gave times = 0:50, it would return 51 results at 0,1,2, ... 50. Other wise I have to supply a sequence such as , times = seq(0,50,0.1).

I have a function which changes rapidly at the beginning and then much more gradually. In MATLAB the output results reflects this with 82 time steps returned in the results, of which 49 are between the time step 0 and 1.

I want to know if there is a way to get R to return the results in the same manner as MATLAB, so without me having the pre-specify the time points I wish the results returned at.

score 2 · Accepted Answer · answered Jun 20 '16 at 09:57

The R-Package deSolve describes the used parameter times in this way:

time sequence for which output is wanted;

The document Dennis liked to has another important sentence:

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.

A simple example is the following one, the Lorenz equations, mentioned in the book about the package 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

You can see the differences in the results at t = 100. This comes from the different defined times.

Regards,
J_F

Is there a way to determine the "best" time interval for the times argument? — Yan, Aug 02 '19 at 17:45

score 1 · Answer 2 · answered Jun 20 '16 at 09:09

Upon reading this document on deSolve it states:

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

As such it would seem that you should indeed use times = seq(0,50,0.1) as input. If you only want to show 'interesting' points in a graph, I suppose you could write a small function to post-process the output of the actual solver output.

ODE Times Matlab vs R

2 Answers2

Linked