[R-sig-dyn-mod] Time steps of integration

I didn't mean to imply that every step it tried would be kept, only the ones that were considered acceptable. I haven't found any consistent terminology, but in Radhakrishnan and Hindmarsh they seem to use "steps" for the acceptable solution, and "iterations" for the repeated attempts to find a solution at each step. Following Daniel's example, modifying the code as such seems to keep steps but not iterations. It is, however, not entirely clear that it actually keeps steps and only steps and there are other problems with such an implementation.

modelfunc = function(curr_time, state, parms) {
	dy = cos(curr_time)
	if( curr_time > gl.t[length(gl.t)] ) {
		gl.t <<- c( gl.t, curr_time)
		gl.s <<- c( gl.s, state) 
	} else {
		gl.t[length(gl.t)] <<- curr_time
		gl.s[length(gl.s)] <<- state
	}

	return(list(dy))
}

times=c(0, 8*pi)  # Note that this only contains the start and end times.
yinitial = 0
names(yinitial) = 'y'

gl.t <- 0
gl.s <- 0

aliasedodeout = lsodes(yinitial, times, modelfunc)

x = seq(0, 8*pi, by=pi/16)

plot(sin(x) ~ x, type='l')
points(gl.s ~ gl.t, col='darkgreen')
abline(v=seq(0, 8*pi, by=pi), lty=2)

with(as.data.frame(aliasedodeout), {
       points(y ~ time, type='b', col='blue')   
})

attributes(aliasedodeout)$istate[2]
attributes(aliasedodeout)$istate[3]
length(gl.t)  # This is 1 longer than the number of steps, but considerably shorter than the number of function evaluations, so it seems to contain the starting value and each successful step of integration.

It appears that gl.t and gl.s contain times and states at steps that the solver thought were appropriate before moving to the next time step, but they do not contain the iterations required at each step. The plot shows that gl.t and gl.s give an excellent representation of the function. Note that I only gave lsodes the start and end times of the interval, and I did not have to think at all about whether the times I specified were appropriate! Also note that the blue line, the result lsodes actually gave me, is worthless, even though it actually calculated a great solution.

In terms of efficiency, simply returning the successful time steps and states could be more efficient in some cases. If a user naively gave a highly dense set of times, the solver makes needless calculations and stores unnecessary data.

times = seq(0, 8*pi, by=pi/64)
yinitial = 0
names(yinitial) = 'y'

gl.t <- 0
gl.s <- 0

noaliasedodeout = lsodes(yinitial, times, modelfunc)

diagnostics.deSolve(noaliasedodeout)

There are more than twice as many function evaluations and values to store. In a situation where the user only cares that the solution is within the error tolerances and does not care about the specific times in the solution (I suspect this common), the first solution is more efficient computationally and spatially.

In Radhakrishnan and Hindmarsh (1993) they make the distinction between xi_out[j], the values of the independent variable at which the user wants a solution, and xi[n], the values of the independent variable that the solver arrived at after multiple iterations at each step. Of course if a user wants specific times, then xi_out[j] is important, but one should always be interested in xi[n] because throwing them out potentially removes information in quickly varying portions of the function. Although I like what Daniel has done, it is not a very elegant solution. I'm not asking that anything be taken away from what currently exists, only to provide the option to keep the times and states at the steps (not iterations) of the solver without having to resort to global variables or other chicanery. It would seem simple enough to return a solution that contained times and states at both xi_out[j] and xi[n]. Perhaps it could be made available to the user by adding an option such as k!
eep_solver_steps=TRUE?

I hope I do not sound like I dislike this package. I think it is exceptionally good, but this would be quite a welcome addition that would solve Lena's problem as well as some others.

Best wishes,
Justin