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

[Karline Soetaert]

Dan,

What you see is the numerical creation of the Jacobian, by perturbation.

 For an implicit method like the "bdf", at each time step each of the state variables is perturbed with a small value (which relates to atoll/rtol), and the derivatives calculated. 
This is to estimate the Jacobian (d(dCi/dt) / dCj) which is used to advance the simulation to the next step.

When you use "lsoda" the model is solved explicitly, hence no jacobian is needed.

If you do: plot(ode(y=C,times=seq(0,50,by=10),func=Model,parms=p,method="bdf"))

Then you will find that the result is ok, even when using the bdf method.



Karline
-----Original Message-----
From: r-sig-dynamic-models-bounces at r-project.org [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Daniel Reed
Sent: dinsdag 24 april 2012 16:30
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] Solving ODEs using the BDF method

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!
_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models