[R-sig-dyn-mod] hmax and its role in ode

[Thomas Petzoldt]

Andras,

your rlnorm() distribution may produce relatively high values of Ka, and
processes with high turnover rates are well-known to produce numerical
problems in general. The reason is the limited accuracy of "ordinary"
software:

The smallest absolute value of a machine can be obtained with:

.Machine$double.xmin

and the minimum (relative) accuracy with:

.Machine$double.eps

There are several methods to handle such issues. All these methods
require special care and good knowledge of the underlying processes and
the numerical limitations of computers.

One method is to rescale the problem to another time unit (e.g. from
daily to hourly rates), resp. to decrease hmax manually.

As another option you may try to transform your equations to log-scale.
The latter was discussed some time ago on this list, see also the
example below.


Thomas P.




library(deSolve)

derivs     <- function(t, y, p) list(y = -p * y)

derivs_log <- function(t, y, p) list(y = - p)


p <- 2  # change this to 0.1, 2, 4, 5, 10


## accuracy is limited by .Machine$double.eps
out <- ode(y = c(y = 600), times = 1:150, atol=1e-15, rtol=1e-15,
   func = derivs, parms = p, hmax = 1e-2)

## log transformed is more robust, but limited by .Machine$double.xmin
out_log <- ode(y = c(y = log(600)), times = 1:150,
   func = derivs_log, parms = p)

plot(out, mfrow = c(2, 1), lwd=3)
lines(out_log[,1], exp(out_log[,"y"]), col="green",
   lwd=2, lty="dashed")

plot(out_log[,1], exp(out_log[,"y"]), log="y", col="green",
   type="l", lwd=3, main="log(y)")
lines(out[out[,"y"] > 0, ], lwd=2, lty="dashed")

# compare results and check for negative values
summary(out[,"y"] - exp(out_log[,"y"]))
min(out[,"y"])
min(exp(out_log[,"y"]))

## accuracy of the machine
.Machine$double.xmin
.Machine$double.eps