[R-sig-dyn-mod] ode integration troubles

[Thomas Petzoldt]

Hi Romain,

I've made quick test, and you are right, the model runs only to a
certain point (time step 63) and then produced only NAs. This behavior
is typical for models where some functions are not protected against
"wrong" values, e.g. if values become negative.

(1) In my opinion, using the "safeguard-method" of setting values to
max(X, 0) is the wrong way. This approach potentially violates mass
balance and it makes the model difficult for the solvers.

You should always keep in mind, that the model function returns the 
derivatives, not the new states!

Instead of setting hard limits for the *state variables*, I would
recommend to  make sure that a quantity that becomes negative is
safeguarded by the functional form of its *turnover rates*. An example
how this can be done is given below. The half sat constant (here kloss)
is uncritical as long as it is small (compared to the state variable,
but not too small to avoid numerical problems.

(2) Another point, that I found is that your state variables have very
different orders of magnitude. In such cases, it may be necessary to set
the tolerance parameters (atol and rtol) to individual values for each
state variable (as vectors). The default for both is 1e-6.

Often, relative tolerance (rtol) can remain the same, whereas absolute
tolerance has (atol) has to respect order of magnitude, e.g. 1e-6 if a
value is around 1 (e.g. between 0.01..100) and 1e-12 if the value is
around 1e-6.

If I look at plot(out) of your model, then potentially dangerous
behavior can be seen for Lysis_ZL, PO4_uptake_Ulva, PO4_uptake_graci,
f_QP_Zost ZGR, f_LP, f_SN_zost, LABS_PO4, FxNO3, ... and maybe also for
the nutrients (PO4, NH4, NO3).


(3) A way to debug your model is to switch "candidate processes" off
resp. to set them to fixed values. If the model runs through, enable
them step by step.

Best regards,

Thomas




## without safeguard
modelA <- function (time, y, parms) {
   with(as.list(c(y, parms)), {
     f <- S/(kS + S)
     dX <- r * f * X
     dS   <- - r * f * X - loss
     list(c(dX, dS))
   })
}



y      <- c(X = 1, S = 1)
parms  <- c(r = 0.1, kS = 0.5, loss=0.01)

outA1 <- ode(y, 0:100, modelA, parms)
plot(outA1) # looks o.k., but is already wrong !!!

outA2 <- ode(y, 0:365, modelA, parms)
plot(outA2) # obviously wrong


## with "Monod-type" safeguard
modelB <- function (time, y, parms) {
   with(as.list(c(y, parms)), {
     f <- S/(kS + S)
     loss <- maxloss * S / (kloss + S) # <-- safeguard
     dX <- r * f * X
     dS   <- - r * f * X - loss
     list(c(dX, dS))
   })
}

parms  <- c(r = 0.1, kS = 0.5, maxloss=0.01, kloss=1e-4)

outB <- ode(y, 0:365, modelB, parms)

## compare A and B
plot(outA1, outB)