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

Mon Nov 18 09:30:42 CET 2013

Hi,

Just a few remarks.

First, none of the solvers in deSolve return the time step nor do they store the values of state variable values at all time steps during the simulation. 
For the odepack solvers, there is a history array that keeps the solution and derivatives of the most recent times, and which is stored in a common block.  This is a rather small matrix, as it keeps only the last values. The solvers from the radau family store a vector with polynomial coefficients. 
Thus, returning the values of state variables at all retained time steps would require serious hacking in the underlying code. A complicating factor is that we do not know in advance how many time steps will be taken.

Secondly, the value at the requested output times is not found by simple linear interpolation at the end (as suggested by Justin), but higher-order interpolation during the integration, usually with variable order, the order depending on the dynamics. Thus, allowing to let the user to linearly interpolate afterwards would give much less accurate solutions. 

As for the more sophisticated advection schemes (Daniels request), there are a few in the ReacTran package, and they are to be used with constant time steps. For integration methods that adapt the time step, the actual size of the time step is only known after the step has been decided upon by the integration solver. Thus, for using a advection scheme that depends on the time step, one would need to do some kind of iteration.

Karline

-----Original Message-----
From: r-sig-dynamic-models-bounces at r-project.org [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Daniel Reed
Sent: vrijdag 15 november 2013 22:41
To: Special Interest Group for Dynamic Simulation Models in R
Subject: Re: [R-sig-dyn-mod] Time steps of integration

Hi all:

Such functionality would also be useful when implementing numerical methods (e.g., advection schemes) that require knowledge of the current time step. If the ODE solvers had the option of keeping tabs on time steps, this would certainly help and would really add to what is already a stellar package. 

Cheers,
Daniel

_______________________________
Daniel C. Reed, PhD.
Postdoctoral Fellow,
Dept. of Earth & Environmental Sciences, University of Michigan, Ann Arbor, MI, USA.
email: reeddc at umich.edu
web: www.danielreed.org

On Nov 15, 2013, at 3:39 PM, "McGrath, Justin M" <jmcgrath at illinois.edu> wrote:

> 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
> 
> ________________________________________
> From: r-sig-dynamic-models-bounces at r-project.org 
> [r-sig-dynamic-models-bounces at r-project.org] on behalf of Thomas 
> Petzoldt [Thomas.Petzoldt at tu-dresden.de]
> Sent: Friday, November 15, 2013 11:08 AM
> To: Special Interest Group for Dynamic Simulation Models in R
> Subject: Re: [R-sig-dyn-mod] Time steps of integration
> 
> Hi,
> 
> the solver interpolates *during* the integration. Saving all 
> intermediate steps would be much too inefficient and memory consuming.
> 
> Sometimes, considerable work is needed to identify a good automatic 
> stepping and sometimes, steps are discarded and the solver goes a 
> little bit back.
> 
> In case of events, time steps are calculated this way, that events are 
> matched. If no events are present, integration is allowed to overrun 
> the external steps.
> 
> If you have external forcings, step size should be limited with hmax.
> 
> 
> Thomas
> 
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> 
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models