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

[Thomas Petzoldt]

On 04.02.2016 10:35, Romain PETE wrote:
> 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)

What I mean is that the problem results possibly from an incomplete
model specification. This can be due to a programming error, missing
processes or of inadequate functional descriptions of a process.

Let's assume you have a state variable, e.g. a nutrient that needs
to be > 0, than all loss processes have to be formulated in a way
that no loss can occur if the state goes to zero.

The introduction of such a "soft safeguard" (or better negative
feedback) is only necessary for processes, missing adequate feedback
mechanisms. A Monod term with a low half saturation constant means
essentially that a rate is mostly constant, except shortly before a
state (e.g. nutrient) goes to zero.

Constructions like rate = max(rate, eps) with eps > 0 (e.g. 1e-6 or
whatever) should always be a very last resort measure. Such hard limits
can be a problem for the numerical solvers, and lead often to mass
balance errors. Limiting a state variable is even worse.

> 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)

O.k., "candidate" was just an explanation and no scientific term. It
means to use biological knowledge in such cases and to think about which
curves make no sense from a scientific viewpoint. This should then allow
to identify and correct the process causing implausible behavior.

My second mail shows you how to track this down, but you can also
consider to output additional intermediate values or just use print
statements in your function.

Good luck!

Thomas