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

[Setzer.Woodrow at epamail.epa.gov]

Hi, all.
lsoda can certainly not handle complex parameters.  You can try (as Hank
suggested) limiting hmax.  You can also crank up relative and absolute
precision by specifying smaller values of  rtol and atol.  I've seen
similar problems in which the state variable becomes negative, with very
small absolute value, when theoretically, the system has a non-negative
solution.  This is certainly due to imprecision in the numerical
solution.  Have you tried including an analytic jacobian?  That could
improve the numeric properties of the solution.

Woody

R. Woodrow Setzer, Ph. D.
National Center for Computational Toxicology
US Environmental Protection Agency
Mail Drop B205-01/US EPA/RTP, NC 27711
Ph: (919) 541-0128    Fax: (919) 541-1194


                                                                        
             "Martin Henry H.                                           
             Stevens"                                                   
             <HStevens at muohio                                        To 
             .edu>                    Spencer Graves                    
                                      <spencer.graves at pdf.com>          
             06/11/2007 01:02                                        cc 
             PM                       Jeremy Goldhaber-Fiebert          
                                      <JGOLDHAB at hsph.harvard.edu>,      
                                      R-Help                            
                                      <r-help at stat.math.ethz.ch>,       
                                      Woodrow Setzer/RTP/USEPA/US at EPA   
                                                                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"