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

[Romain PETE]

Dear Thomas,

Many thanks for this very quick and informative reply!
I will investigate this safeguard option on turn over rates: if I 
understand correctly, any loss in my rates of change may introduce 
negative values and this is where I need safeguards? (mostly for those 
cited below)

One last question about point3, what you're calling "candidate 
processes" are these processes identified with potentially dangerous 
behaviour? (e.g. decreasing rapidly toward zero or negative values)

Many thanks again indeed

Cheers
romain
Le 03/02/2016 19:28, Thomas Petzoldt a écrit :
> 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)
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

-- 
Romain Pete
UMR MARBEC
Ifremer - Laboratoire LER/LR
Avenue Jean Monnet
CS 30171 - 34203 Sète cedex France
Tél : (+33)4 99 57 32 93 - Fax : 04 99 57 32 96