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



______________________________________________
R-help@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.

Reply via email to