[R-sig-dyn-mod] Problem with simecol and ode - crash from R

[Thomas Petzoldt]

Soetaert, Karline wrote:

[...]

 > A more appropriate model would be:
 >
 >
 > NoCrash <- function (time, N, parms) {
 >             with(as.list(parms),
 >                  {
 >                    dn <-  r * N * (1 - (N/K)) - h*(N>0)
 >                    if (! is.numeric(dn)) dn <- 0
 >                    list(dn)
 >                  }
 >                  )
 >           }
 >
 > require(deSolve)
 > parms  = c(r = 0.1, K = 10000, h=200)
 > init   = c(N = 2000)

Well, this works in practice, but allow me to two additional points:

First, the line "if (! is.numeric(dn)) dn <- 0" is redundant, probably a 
residual of testing?

Second, instead of a hard limit (that is of course a good test case for 
the solvers), I usually prefer a more smooth limit, e.g:

NoCrash <- function (time, N, parms) {
              with(as.list(parms), {
                    dn <-  r * N * (1 - (N/K)) - h * N / (N + hlim)
                    list(dn)
              })
            }

parms  <- c(r = 0.1, K = 10000, h = 200, hlim = 1)
init   <- c(N = 2000)

out <- ode(func=NoCrash, y = init, parms = parms,
   times = seq(0,50,.01), method = "lsoda")

plot(out, type = "l")

You can of course try different values of hlim, e.g. hlim=100.

Thomas P.