1

I am writing a numerical model in R, for an ecological system, and solving it using "lsoda" from package deSolve. My model has 14 state variables. I define the model, set it up fine, and give time duration according to this:

nyears<-60

ndays<-nyears*365+1

times<-seq(0,nyears*365,by=1)

Rates of change of state variables (e.g. the rate of change of variable "A1" is "dA1")are calculated according to existing values for state variables (at time=t) and a set of parameters. Simplified example:

dA1<-Tf*A1*(ImaxA*p_sub)

Where Tf, ImaxA and p_sub are parameters, and A1 is my state variable at time=t.

When I solve the model, I use the lsoda solver like this:

out<-as.data.frame(lsoda(start,times,model,parms))

Sometimes (depending on my parameter combinations), the model run completes over the entire duration I have specified, however sometimes it stops short of the mark (still giving me output up until the solver "crashes"). When it "crashes", this message is displayed:

DLSODA-  At current T (=R1), MXSTEP (=I1) steps
      taken on this call before reaching TOUT
In above message, I1 = 5000

In above message, R1 = 11535.5

Warning messages:
1: In lsoda(start, times, model, parms) :
  an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxsteps
2: In lsoda(start, times, model, parms) :
  Returning early. Results are accurate, as far as they go

It commonly appears when one of the state variables is getting exponentially bigger, or is tending very near to zero, however sometimes it crashes when seemingly not much change is happening. I may be wrong, but is it due to the rate of change of state-variables becoming too large? If so, why might it also "crash" when there is not a fast rate of change? Is there a way that I can make the solver complete its task with the specified parameter values, maybe with a more relaxed tolerance for error?

JB8365
  • 21
  • 1
  • 4
  • no real help, but seems to be a question for stats.stackexchange.com what do you mean with _build in_? Is it `lsoda` from package 'deSolve' ? – ckluss Jan 04 '15 at 15:45
  • Some differential equations simply do not have a stable solution. What level of mathematical sophistication do you possess? (At the moment it's not possible to discuss this with any specificity because there is so much missing information, so I'm voting to close on SO for lack of clarity. As pointed out above, this question might be welcomed, when augmented by further details, at CrossValidated.com) – IRTFM Jan 04 '15 at 18:12
  • it's hard to know without more detail. One thing that can be done to help is to reframe some of the equations in terms of per capita rates, switching from `dX/dt = X*(Y+Z)` to `d(log X)/dt = Y+Z` (you then have to use `exp(log(X))` in any other contexts in the gradient computation) – Ben Bolker Jan 04 '15 at 18:39

1 Answers1

1

Thank you all for your contributions. I looked at some of the rates, and at the point of crashing, the model was switching between two metabolic states - and the fast rate of this binary switch caused the solver to stop - rejecting the solution because the rate of change was too large. I have fixed my model by introducing a gradual switch between states (with a logistic curve) instead of this binary switch. I aknowledge that I didn;t give enough info in the original question, so thanks for the help you offered!

JB8365
  • 21
  • 1
  • 4