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

[Thomas Petzoldt]

Hello,

the cause for the misbehavior lies in the, compared to the state
variables, much smaller half saturation constant of the Monod function
K=1e-6. This is (at least to me) well known, so I always remind users to
use not too small half sat values, especially if they are used as
safeguards preventing state variables to become negative.

A little bit more complex (ecological) models crash in such
circumstances, while Radovans example skips the singularity. Reducing
the tolerances would not solve it, as one may then set a K = 1e-9 or
whatever. That's why we always remind users to be careful and aware of
such pitfalls.

I could not resist and tested the model with the system dynamics tool
"Simile" (www.simulistics.org) -- it shows the same behavior, so I
wonder if and how other "different software" handles this.

Any ideas?

Thomas


On 29.04.2016 16:02, Radovan Omorjan wrote:
> 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
>
> =======
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models


-- 
Dr. Thomas Petzoldt
Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany

E-Mail: thomas.petzoldt at tu-dresden.de
http://tu-dresden.de/Members/thomas.petzoldt

-- limnology and ecological modelling --