Re: [R] Using odesolve to produce non-negative solutions

From: Thomas Petzoldt <thpe_at_simecol.de>
Date: Thu, 14 Jun 2007 11:03:58 +0200

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 ...]



R-help_at_stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Thu 14 Jun 2007 - 09:08:17 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 14 Jun 2007 - 10:32:00 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.