[R-sig-dyn-mod] Solving ODEs using the BDF method

[Daniel Reed]

Hi there:

I have a question regarding the BDF method of solving ODEs, as implemented in the deSolve package. Included below is a trivial example of a model:

library("deSolve")

Model<-function(t,C,p){
        print(C)
        return(list(rep(0,p)))
}

p<-100
C<-rep(0,p)

ode(y=C,times=seq(0,50,by=10),func=Model,parms=p,method="lsoda")

This runs as expected and the model solution is zero throughout over the given time domain. Nevertheless, when I run the same code using BDF as the method instead of LSODA I get a curious result: the number 1e-6 is introduced to the solution. The first three function calls are fine (zeroes throughout C), but on the fourth function call 1e-6 appears as the first element of C. At each of the following function calls 1e-6 moves along to the adjacent element. The erroneous value appears to be linked to the parameter atol. That is, if I set atol as 2, then the value is equal to 2 (same behaviour as before). Clearly, I'm doing something wrong when using BDF, but I cannot work out what it is. Could someone enlighten me?

Many thanks in advance!