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

From: Ravi Varadhan <rvaradhan_at_jhmi.edu>
Date: Mon, 11 Jun 2007 13:23:19 -0400

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.

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 Mon 11 Jun 2007 - 19:52:38 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 Mon 11 Jun 2007 - 20:33:41 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.