[R] Using odesolve to produce non-negative solutions

[Thomas Petzoldt]

Dear Jeremy,

a few notes about your model: The equations of your derivatives are 
designed in a way that can lead to negative state variables with certain 
parameter combinations. In order to avoid this, you are using "if 
constructions" which are intended to correct this. This method is 
however (as far as I have learned from theory and own experience ;-) a 
bad idea and should be strictly avoided.

This trick may violate assumptions of the ODE solvers and in many cases 
also mass-balances (or similar).

Instead of this, processes should be modeled in a way that avoids 
"crossing zero", e.g. exponential decay (x, k > 0):

dx = - k * x  (1)

and not a linear decay like

dx = -k       (2)

which by its nature can lead to negative state values.

Case (1) can be managed almost perfectly by lsoda with his automatic 
internal time step algorithm. hmax is intended for non-autonomous models 
to ensure that external signals are not skipped by the automatism, which 
may be appropriate in your case because p seems to contain time 
dependent functions.

As far as I can see, your equations belong to type (1) and should not 
need any of the if and for constructs, as long as your parameters (which 
are not given in your post) do have the correct sign and range (for 
example, vax <= 1, death >= 0  etc.).

If you perform optimization, you must ensure that parameters stay in the 
plausible range is met using transformations of the parameters or 
constraints in the optimization procedure.

Thomas

PS: another question, what is the purpose of the state variable N?
I guess it can be derived from the other states.


Jeremy Goldhaber-Fiebert wrote:
> Hello,
> 
> I am using odesolve to simulate a group of people moving through time
> and transmitting infections to one another.
> 
> In Matlab, there is a NonNegative option which tells the Matlab
> solver to keep the vector elements of the ODE solution non-negative
> at all times. What is the right way to do this in R?
> 
> Thanks, Jeremy

[... code deleted, see original post ...]