Adaptive time step for Runge-Kutta method ode45 in deSolve

Question

I want to use the explicit Runge-Kutta method ode45 (alias rk45dp7) from the deSolve R package in order to solve an ODE problem with variable step size.

According to the deSolve documentation, it is possible to use adaptive or variable time steps for the rk solver function with the ode45 method instead of equidistant time steps but I'm at loss how to do this.

The rk funcion is called like this:

rk(y, times, func, parms, rtol = 1e-6, atol = 1e-6, verbose = FALSE, tcrit = NULL,
hmin = 0, hmax = NULL, hini = hmax, ynames = TRUE, method = rkMethod("rk45dp7", ... ), 
maxsteps = 5000, dllname = NULL, initfunc = dllname, initpar = parms, rpar = NULL, 
ipar = NULL, nout = 0, outnames = NULL, forcings = NULL, initforc = NULL, fcontrol = 
NULL, events = NULL, ...)

with times being the vector of times at which explicit estimates for y are desired.

For equdistant time steps with a distance of 0.01 I can write times as

times <- seq(0, 100, 0.01)

Supposing I want to solve the equation for the interval from 0 to 100 how do I define times without giving a step size?

Any help would be greatly appreciated.

Answer 1

There are two issues here. First, if you want to specify a vector of times with multiple increments, use this (for example):

times <- c(seq(0,0.9,0.1),seq(1,10,1))
times
#  [1]  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0 10.0

Here, we have [0,1] by 0.1, and [1,10] by 1.

But you don't actually have to do this: the parameter times= tells rk(...) for which times to report results. An adaptive algorithm will internally adjust the time increment to yield accurate results at the times specified in the argument. So for an adaptive algorithm, e.g., method="rk45dp7" you don't have to do anything. For a non-adaptive algorithm, for example method="euler", the time increment used by the algorithm is indeed the increment specified in times=. You can see the effect of this in this simple example, which integrates the Van der Pol oscillator.

y.prime <- function(t,y.vector,b) {    # Van der Pol oscillator
  x <- y.vector[1]
  y <- y.vector[2]
  x.prime <- y
  y.prime <- b*y*(1-x)^2 - x
  return(list(c(x=x.prime,y=y.prime)))
}
h  <- .001                  # time increment
t  <-  seq(0,10,h)          # times to report results
y0 <- c(0.01,0.01)          # initial conditions
euler   <- rk(y0, t,func=y.prime,parms=1,method="euler")
rk45dp7 <- rk(y0, t,func=y.prime,parms=1, method="rk45dp7")
# plot x vs. y
par(mfrow=c(1,2))
plot(euler[,2],euler[,3], type="l",xlab="X",ylab="Y",main=paste("Euler: h =",format(h,digits=3)))
plot(rk45dp7[,2],rk45dp7[,3], type="l",xlab="X",ylab="Y",main=paste("RK45DP7: h =",format(h,digits=3)))

Below is a comparison of the results for several values of h. Note that with method="rk45dp7" the results are stable for h < 0.5. This is because rk45dp7 is internally adjusting the time increment as needed. For method="euler" the results do not match rk45dp7 until h~0.01.