[R-sig-dyn-mod] deSolve::ode workaround for solutions with NaN

[Thomas Petzoldt]

On 30.06.2010 11:32, Maik Renner wrote:
> Dear all,
>
>
>
> thank you a lot for your help and suggestions.
>
> I tried Stephens approach with some success and it was the solution I
> was looking for.
>
> I still have to play around with the values I set for dy.
>
> Because setting it to zero, results in zeros for the remaining
> simulation. But this is problem specific.
>
> Maik


Dear Mike,

yes, you are right, this is strictly problem specific. Crossing zero is 
a problem of your model definition and its parametrization, and *not* an 
instability of the solver.

The solution of Stephen Ellner works as long as mass balances are not 
violated (e.g. if other states do not 'profit' from the negative 
derivative before it is zeroed).

If the scientific question allows, then it is better to adapt the model 
itself than tuning the tools. So I would prefer something like a 
Michaelis-Menten function instead of a constant loss, e.g. something 
like the following:


library(deSolve)

func1 <- function(t, y, p) {
   with(as.list(p), {
    dy <-  -a
    list(c(dy), c(a))
   })
}

func2 <- function(t, y, p) {
   with(as.list(p), {
    a  <- y * b / (k + y)
    dy <-  -a
    list(c(dy), c(a))
   })
}


times <- seq(0, 10, .1)
y <- 5

o1 <- ode(y, times, func1, c(a = 1))

o2 <- ode(y, times, func2, c(b = 1, k = 1e-6), atol=1e-8)

par(mfrow=c(1,2))

matplot(o1[,1], o1[,-1], type="l", ylim=c(-1,5))
abline(h=0, lty="dotted")

matplot(o2[,1], o2[,-1], type="l", ylim=c(-1,5))
abline(h=0, lty="dotted")


The half saturation parameter k is uncritical but it should
not be too low because of limited numerical precision. If you set it to 
a 'very' small number (e.g. 1e-6), you need to adjust atol accordingly 
(e.g. 1e-8). Instead of a MM you can also use a sigmoid function.



Thomas Petzoldt