[R-sig-dyn-mod] stiff ode - accuracy?

[Radovan Omorjan]

Hello,

Here is one of my test examples for stiff ode solvers (see the code 
below). I tried many lsoda based solvers in different software and they 
usually can solve this particular problem with their default parameter 
values (accuracy etc.). By using ode() from deSolve, it seems that I 
have to reduce the values of atol and rtol in order to get the expected 
profiles. There is no method using the default accuracy values I tried 
which will give me the profile I expect.

The problem with these equations become problematic when times riches 
around 16.33. After that time the value of B should be constant and for 
S should be zero.

QUESTION: I just wanted to send this example to the mailing list and 
hope there might be any comment about it. In the first place, I hope 
that the authors of deSolve might consider this example regarding the 
default values of atol and rtol.

Best Regards,
Radovan

==========

library(deSolve)

par <- c(k = 0.3, K = 1e-6)

yini <- c(B = 0.05, S = 5)

biomass <- function (t, y, parms) {
   with(as.list(c(y,parms)), {
     dBdt <- k*B*S/(K + S)
     dSdt <- -0.75*k*B*S/(K + S)
     list(c(dBdt, dSdt))
   })
}
times <- seq(from = 0, to = 20, by = 0.01)

#NO SUCCESS
out<- ode(y = yini,
           times = times,
           func = biomass,
           parms = par,
           method = "lsode") #we can use any available method
plot(out, lwd = 2)
diagnostics(out)

#If we just lower the tolerance, the picture is different
# the profile for B should remain constand and for S should be zero

out<- ode(y = yini,
           times = times,
           func = biomass,
           parms = par,
           rtol = 1e-7,
           atol = 1e-7)
plot(out, lwd = 2)
diagnostics(out)
##End

=======