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

[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')	
})

>From the plot, the solver seems to have reached the proper solution in each case, but for the blue line, the fact that I've asked for fewer times than required to represent this function properly has resulted in aliasing. I.e., the solver did everything right in each case, the problem lies in the fact that I didn't get the solution the solver arrived at, I got a subset of it.

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
________________________________________
From: r-sig-dynamic-models-bounces at r-project.org [r-sig-dynamic-models-bounces at r-project.org] on behalf of Daniel Reed [reeddc at umich.edu]
Sent: Tuesday, November 05, 2013 6:20 AM
To: Special Interest Group for Dynamic Simulation Models in R
Subject: Re: [R-sig-dyn-mod] Time steps of integration

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:

> Dear all,
>
> I face some trouble trying to figure out the real step size of the integration, when using the ode solver (radau, ode23... etc). I get a summary with verbose telling me how many time steps it used in total, but I would rather need a vector with all time steps.
> Is there also a possibility to get the time needed for the integration, although I think I can handle this with the help of proc.time or system.time.
>
> Thanks a lot in advance,
> Lena
>
>       [[alternative HTML version deleted]]
>
> _______________________________________________
> 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