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

Daniel Reed reeddc at umich.edu
Tue Apr 24 16:59:33 CEST 2012


That makes a lot of sense now. Thanks for the insight, Karline!


On 24 Apr 2012, at 10:47, Karline Soetaert wrote:

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



More information about the R-sig-dynamic-models mailing list