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

[Ravi Varadhan]

Hi Jeremy,

A smaller step size may or may not help.  If the issue is simply truncation
error, that is the error involved in discretizing the differential
equations, then a smaller step size would help.  If, however, the true
solution to the differential equation is negative, for some t, then the
numerical solution should also be negative.  If the negative solution does
not make sense, then the system of equation needs to be examined to see when
and why negative solutions arise.  Perhaps, I am just making this up - there
needs to be a "barrier function" that slows down the trajectory as it
approaches zero from its initial value. It is also possible that only
certain regions of the parameter space are allowed in the sense that only
there the solution is feasible for all t.  So, in your example, the
parameters might not be realistic.  In short, if you are sure that the
numerical solution is accurate, then you need to go back to your system of
equations and analyze them carefully.


Ravi.


----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 

----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Martin Henry H.
Stevens
Sent: Monday, June 11, 2007 1:03 PM
To: Spencer Graves
Cc: Jeremy Goldhaber-Fiebert; R-Help; Setzer.Woodrow at epamail.epa.gov
Subject: Re: [R] Fwd: Using odesolve to produce non-negative solutions

Hi Spencer,
I have copied Woody Setzer. I have no idea whether lsoda can estimate  
parameters that could take imaginary values.
Hank
On Jun 11, 2007, at 12:52 PM, Spencer Graves wrote:

> <in line>
>
> Martin Henry H. Stevens wrote:
>> Hi Jeremy,
>> First, setting hmax to a small number could prevent a large step, if
>> you think that is a problem. Second, however, I don't see how you can
>> get a negative population size when using the log trick.
> SG:  Can lsoda estimate complex or imaginary parameters?
Hmm. I have no idea.
>
>> I would think that that would prevent completely any negative values
>> of N (i.e. e^-100000 > 0). Can you explain? or do you want to a void
>> that trick? The only other solver I know of is rk4 and it is not
>> recommended.
>> Hank
>> On Jun 11, 2007, at 11:46 AM, Jeremy Goldhaber-Fiebert wrote:
>>
>>> Hi Spencer,
>>>
>>> Thank you for your response. I also did not see anything on the  
>>> lsoda
>>> help page which is the reason that I wrote to the list.
>>>
>>>> From your response, I am not sure if I asked my question clearly.
>>>
>>> I am modeling a group of people (in a variety of health states)
>>> moving through time (and getting infected with an infectious
>>> disease). This means that the count of the number of people in each
>>> state should be positive at all times.
>>>
>>> What appears to happen is that lsoda asks for a derivative at a  
>>> given
>>> point in time t and then adjusts the state of the population.
>>> However, perhaps due to numerical instability, it occasionally lower
>>> the population count below 0 for one of the health states (perhaps
>>> because it's step size is too big or something).
>>>
>>> I have tried both the logarithm trick
> <snip>
>
> ______________________________________________
> 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.



Dr. Hank Stevens, Assistant Professor
338 Pearson Hall
Botany Department
Miami University
Oxford, OH 45056

Office: (513) 529-4206
Lab: (513) 529-4262
FAX: (513) 529-4243
http://www.cas.muohio.edu/~stevenmh/
http://www.muohio.edu/ecology/
http://www.muohio.edu/botany/

"E Pluribus Unum"

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