Skip to content

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

9 messages · Lena Appel, Daniel Reed, McGrath, Justin M +2 more

Original
Daniel Reed
# 
Hi Lena:

That?s a good question. I?m sure one of the authors of deSolve can give you a better solution, but here are my thoughts based on my limited understanding of how the ODE integrators work. One approach is to define a global array that stores step times, say, gl.t. So, before the model solution you could define your start point

#Start time is zero
gl.t <- 0

Then, in your model function, have some code like this:

if( t > gl.t[length(gl.t)] ) gl.t <<- c( gl.t, t)
else gl.t[length(gl.t)] <<- t

This (untested) piece of code checks to see if the current time used by the ODE solver, t, is later than the previously recorded time. If so, it tags it on to the end of the gl.t array. If not, in the case where ODE solver is trying a smaller time step for the same point in time, the previously recorded value is overwritten. After solving the model, you can determine the time steps using,

timesteps <- diff(gl.t)

As I say, this is an untested idea, but hopefully it helps.

Cheers,
Daniel

PS There is a function called timestep() that I thought was for this very purpose, but it seems to behave differently to how I anticipated. I think this was briefly discussed on this mailing list a while ago.

_______________________________
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 5, 2013, at 4:57 AM, Lena Appel <jumpingfrog at hotmail.com> wrote:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    
10 days later

    

      
      

        


  

        

      McGrath, Justin M
    

    

      
      #
    

  


  

      
I'm chiming in a bit late on this, but why do these functions not simply return the results from the integration? As I understand it, the solver tries multiple time steps at each stage of the integration, and when certain criteria are met, it keeps that time and moves to the next stage. The result is a sequence of unequally spaced time steps. After everything has been solved, the solver takes the time steps you asked for and linearly interpolates the solution to get values at the exact times you requested. There must be some exception to this for events, but I don't know the details.

It would be nice if there were a way to simply return the time steps used in the integration. Then the user could interpolate the data if desired. This would solve Lena's problem, but also one wouldn't have to be concerned that the times they've asked for are appropriate for the solution - the solver should have already done that.

Here's a toy example showing how asking for improperly spaced time points gives a bad solution, even though the solver appears to have worked properly.

modelfunc = function(curr_time, state, parms) {
	dy = cos(curr_time)  # The solution is sin(curr_time).
	return(list(dy))
}

times = seq(0, 8*pi, by=pi/1.8)  # sin(x) has a frequency of 2*pi, by sampling at less than half the frequency, there will be aliasing.
yinitial = 0
names(yinitial) = 'y'

aliasedodeout = lsodes(yinitial, times, modelfunc)

# Plot an accurate depiction of the solution with some guidelines.
x = seq(0, 8*pi, by=pi/32)
plot(sin(x) ~ x, type='l')
abline(v=seq(0, 8*pi, by=pi), lty=2)

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

times = seq(0, 8*pi, by=pi/16)  # Now sample at well less than half the frequency. There won't be aliasing.
noaliasedodeout = lsodes(yinitial, times, modelfunc)

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

            

      
      
      
    

        
Other solvers provide this capability if you pass them only two times, assumed to be the start and end of the period over which you want a solution. It seems that adding this wouldn't really interfere with the current usage since it seems highly unlikely that someone would only want two time points. However, it would seem to make the current usage unnecessary, since one can easily interpolate the result. Here's an example of how it would work.

start_time = 0
end_time = 8*pi
times = c(start_time, end_time)
lsodes(yinitial, times, modelfunc)  # Returns an array with a row for each successful time step during the integration.

hmax would have to be somehow guessed by the solver, or explicitly supplied by the user in this case. Could this be added to deSolve, or do I not understand something about how these solvers work?The example given is pretty ridiculous since the solution is easily derived and it's easy to work out the times needed, but for complex sets of equations I worry that I'm missing detail in quickly changing areas.

Best wishes,
Justin

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Thomas Petzoldt
    

    

      
      #
    

  


  

      
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

  
  


  


      

    

    

      
      

        


  

        

      McGrath, Justin M
    

    

      
      #
    

  


  

      
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 keep_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

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Daniel Reed
    

    

      
      #
    

  


  

      
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:

        

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    
2 days later

    

      
      

        


  

        

      Karline Soetaert
    

    

      
      #
    

  


  

      
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:

        

            

      
      
      
    

        
k!

            

      
      
      
    

        

      
      
    

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

  
  


  


      

    

    

      
      

        


  

        

      McGrath, Justin M
    

    

      
      #
    

  


  

      
Is the conclusion that there is little value in keeping the time steps and state variables at those steps? It seems to me that it would be worth modifying the odepack code, and the way it's written now appears to be an artifact of limited RAM when the algorithm was developed, which is hardly an issue now. Whether interpolation is performed during or following the integration doesn't seem to preclude storing time steps, and not knowing the number of steps is a common and solvable issue with many algorithms.

I hate to beat a dead horse, but it makes more sense to keep the time steps rather than a set of times that has no bearing on the functions being solved. The time steps are almost by definition where the action is happening, and one would only request other times if there were a specific interest in them. Do you not find the first example in my last e-mail, where only the start and end of the period are given, to be not only easier and more intuitive, but to give a better result than a naively chosen set of times? If the answer is that this won't be done because no one cares or there is no one available, then so be it I guess, but then can you confirm that Daniel's suggestion keeps all time steps and only time steps? That way I could store them myself.

Best wishes,
Justin

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Thomas Petzoldt
    

    

      
      #
    

  


  

      
Dear Justin and others,

many thanks for your thoughts and the really nice example, but I don't
think "variable internal steps" would seamlessly fit in deSolve's
philosophy:

odepack provides are a standard set of solvers. They have been already
modified by us to include forcings and events. However, interpolation
was already built in. We are happy with this, and I do not think, that
it would be wise to rewrite the odepack functions.

The Runge-Kutta solvers are a new implementation from the equations, and
partly inspired by Press et al. It would be possible in principle, to
save all intermediate steps as an option, however, I don't assume that
we will do it:

1) Memory: For bigger and esp. stiff problems, saving all intermediate
steps requires more than a bit more of of memory. Karline mentioned that
we do not know in advance how many internal steps are required. This
would require dynamic variables. This is possible, but would further
increase complexity of the code.

2) More important: Some RK solvers (e.g. ode45) use "dense output", a
mechanism, that interpolates directly from the intermediate
Runge-Kutta steps. This is very efficient and an integrated part of the
algorithm and its underlying mathematics, see Karline's comment:

        
On 18.11.2013 09:30, Karline Soetaert wrote:

            

      
      
      
    

        
As said, the cosine example is nice, but it is not very typical and I
don't think that we will implement saving internal steps in deSolve.

However,  if someone wants, one may start to provide such an algorithm
in an own proof of concept package. R is free and open source software.
Contributions are always welcome :)

Thanks,

Thomas