[R] lsoda warning "too much accuracy requested"

[Thomas Petzoldt]

A short addendum, resulting from an off-list discussion:

The reason why Colleen's code failed was raising a negative base to a 
fractional exponent in the third state equation for certain sets of 
parameters, esp. fractional values of beta.

Old versions of odesolve broke down, and recent versions of deSolve (or 
the deprecated odesolve) simply returned NaN for the third state. In 
order to track such problems down, it is helpful to call the system 
separately, e.g.:

 > model(0, start, parms)
[[1]]
         H        BA         N
0.3702103 1.2229436       NaN

and then isolate the problem in the respective state equation.

Thomas P.



Thomas Petzoldt wrote:
> Hi Colleen,
> 
> this error was not uncommon and is usually a sign of a numerically 
> problematic or wrongly implemented model. Please use package deSolve, 
> the successor of odesolve, that is more robust and has also a bunch 
> alternative solvers for difficult cases.
> 
> I tested your code with deSolve (on R 2.8.0, Windows) and it runs 
> without problems.
> 
> Thomas Petzoldt
> 
> 
> BTW: your system worked also with odesolve, so my question: which 
> versions (R, odesolve) and operating system are you using?
> 
> Colleen Carlson wrote:
>> Dear list -
>>
>>  
>>
>> Does anyone have any ideas / comments about why I am receiving the 
>> following
>> warning when I run lsoda:
>>
>>  
>>
>> 1:  lsoda--  at t (=r1), too much accuracy requested in: lsoda(start, 
>> times,
>> model, parms)
>> 2:       for precision of machine..  see tolsf (=r2) in: lsoda(start, 
>> times,
>> model, parms)
>>
>>  
>>
>> I have tried changing both rtol and atol but without success.  I saw the
>> thread in the R-archive of 11 June 2004 but this has not helped me.
>>
>>  
>>
>> I have built the model in stages and the problem only occurs when the
>> exponent beta in the third DE is anything other than 0 or 1.  If beta 
>> = 0 or
>> 1 then the solver gives me perfectly justifiable results.  Just changing
>> beta to 0.9 or similar causes the problem.
>>
>>  
>>
>> I am still new to R so I am unsure if it is my programming or my
>> understanding of the way lsoda works.
>>
>>  
>>
>> Any comments or input would be welcome.
>>
>> Many thanks
>>
>> Colleen
>>
>> ___________
>>
>>  
>>
>> My code is:
>>
>>  
>>
>> library(odesolve)
>>
>> SI <- 80
>>
>> model <- function(t, x, parms) {
>>
>>             H  <- x[1]
>>
>>             BA <- x[2]
>>
>>             N  <- x[3]
>>
>> with(as.list(parms), {
>>
>>             dHdt <- (b/c)*(((a**c)*((H)**(1-c))-H))
>>
>>             dBAdt <- -(BA*b)*(c0+(c1*SI)-log(BA))/(log(1-((H/a)**c)))
>>
>>             dNdt <-  N*alpha*(((log(1-((H/a)**c)))/b)**beta) - (gamma*BA)
>>
>>             list(c(dHdt, dBAdt, dNdt))
>>
>>             })
>>
>>             }
>>
>> times <- seq(0, 40, 1)
>>
>>  
>>
>> parms <- c(a=(SI*1.258621)-1.32759, b=0.1, c=0.4, c0=4.6012, c1=0.013597,
>> alpha=0.0005, beta=0.5,  gamma=0.01)
>>
>> start <- c(H=0.1, BA=0.1, N=600)
>>
>>  
>>
>> out <- as.data.frame(lsoda(start, times, model, parms))
>>
>>
>>     [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help at r-project.org 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.
> 
>