[R-sig-dyn-mod] ReacTran & runsteady

[Daniel Reed]

Hi all:

I was hoping that someone could shed some light on how runsteady works in the ReacTran package. I have a fairly involved model that is producing errors -- specifically crashing due to NaNs in the solution -- so I produced a toy example to see if I could reproduce this behaviour in a simple setting. 

library(ReacTran)

S<-rep(0,1000)

Model<-function(t,S,parms){
	cat(S,"\n")
	return (list(rep(0,length(S))))
}

soln<-steady(y=S,func=Model,parms=c(0),method="runsteady")$y

The above example is trivial and takes an initial array of one thousand elements with values of zero, then returns a dy/dt of zero for all elements so there should be no change. My understanding is that runsteady simply steps the model forward in time until dy/dt is within user-specified tolerances.  In the first two steps, the model behaves as I expect, returning dy/dt=0 and therefore S remains as zero. Thereafter, it produces NaNs. Could someone explain why? I expect it's because I'm trying to peek at solutions before they're ready, during the integration process, but I'd like to understand why I see this behaviour.

Many thanks in advance.