integrate_const with dense output stepper misses final observer call

[demianmnave]

Hi,

In integrate/detail/integrate_const.hpp, the dense output integrator misses the observer call at the final timestep. I fixed the problem by handling the end point interpolation explicitly:

for(;;) {
  Time zm = stepper.current_time();

  // Advance interpolated time up to zm:
  while(start_time <= zm) {
    stepper.calc_state(start_time, start_state);
    obs(start_state, start_time);
    start_time += dt;
  }

  // Advance zm up to end_time:
  if(zm + stepper.current_time_step() < end_time) {
    stepper.do_step( system );
    ++ count;
  }

  // Now handle end_time:
  else break;
}

// Step exactly to end_time:
stepper.initialize(stepper.current_state(),
    stepper.current_time(), end_time - stepper.current_time());
stepper.do_step( system );
++count;

// Step interpolated time as close to end_time as possible without
// crossing it:
while(start_time <= end_time) {
  stepper.calc_state(start_time, start_state);
  obs(start_state, start_time);
  start_time += dt;
}

// Finally, interpolate at end_time if we're far enough away:
if(start_time - end_time < dt/2) {
  stepper.calc_state(end_time, start_state);
  obs(start_state, end_time);
}

// Done:
return count;

I tested this with t=[0,1], and dt=1/4 (no roundoff), dt=.1 (final start_time is less than 1 due to roundoff), and dt=.01 (final start_time is more than 1 due to roundoff). Without the last if(...), and if dt is far from a multiple of the time span, the last two timesteps can be separated by nearly 2*dt. For example, for t=[0,1] and dt=.011, the last two timesteps are separated by a difference of .021. That last if(...) reduces the max difference to dt. Without additional knowledge of the Time type (e.g. via std::numeric_limits()), there doesn't seem to be much more we can do to ensure a reasonable separation between the last two timesteps.

Hopefully you find this code useful.

Cheers, Demian

[headmyshoulder]

Hi Demian, thanks for your feedback. For me it looks like one additional call of the observer after line 121 in integrate/detail/integrate_const.hpp should work:

obs( stepper.current_state() , stepper.current_time() );

The stepper is re-initialized before that line such that the first step will exactly integrate the state to end_time.

On the other hand the function should call the observer in equidistance steps. Imagine for example you record the data inside the observer to use it later in a fft where an equidistant spacing is needed.

Do you need to call the observer at the very end of your integration?

[mariomulansky]

hi demian,

you are pointing out a delicate issue here. The integrate_const routine is supposed to make observer calls after at times t_start + ndt. This means that generally observer calls at the end point are not required except the time span is an exact multiple of dt. The problem now is that due to finite precision we can not tell if it is an exact multiple or not. So in the current implementation the behavior is undefined depending on the rounding behavior. If you want an observer call at your endpoint t, best practice is to use t+dt/2 as endpoint and everything works fine . However, this behavior is obviously confusing and your post convinced that some workaround is needed. So I decided to change the integrate_const function such that it will integrate up to end_time + dt/4. This means the endpoint will be part of the intervall and also if one uses t+dt/2 as endpoint the behavior will still be as expected. However, this can lead to an observer call beyond the originally given end_time which again is confusing. But I assume this case to be less likely so I hope this is a reasonable solution. I will push the fix shortly and I highly appreciate any further comments.

Best, Mario

(_) This actually is exactly how integrate_n_steps works, it calls integrate_const with end_time = start_time + n_dt + dt/2

[demianmnave]

After reviewing the documentation again, I think it's more confusing to force endpoint integration within integrate_const(). The problem with using t1+dt/2 as the endpoint is that t1 is interpolated rather than computed as precisely as possible if the stepper were evaluated exactly at t1 (which is what my version of the function tried to do).

I think the integrate_times() function would actually do the job nicely, albeit with some extra work on the user's part. Maybe a good "fix" would be to implement integrate_n_steps() for Controlled Error and Dense Output steppers in terms of integrate_times()? For example, it could be useful to have an overload of integrate_n_steps() that takes t0, t1, and n_steps as parameters, with dt used to specify the initial timestep.

For my problem (tracing orthogonal projections of curves onto surfaces), I already know the value of the ODE at the endpoint, so having this functionality is not critical. I was actually trying to impose endpoint integration in order to test the accuracy of some of the integration functions for my problem (Bulirsh-Stoer for the current experiment).

Anyway, thanks very much for responding so quickly, and definitely thanks for producing such a useful library.

Cheers, Demian

[mariomulansky]

you are right, my work-around now forces interpolation of the endpoint. even more if the ode had a pole at some t_end+eps my workaround might result in NaN although the user took care not to cross the pole. Your suggestion is right, the root of evil here is the definition of both t_end and dt and disappears when using appropriate integrate_n_steps functions. We will look into this, thanks a lot for your comments.

By the way, can you give us some references on your problem where the endpoints are known? We actually are searching for a way how to test the numerical methods and this sounds as it might be helpful.

Thanks,

Mario

[headmyshoulder]

In my opinion the requirements on integrate_const should be (regardless of the stepper type=

• the observer is called at equidistant steps
• the last step must be at t < t_end + eps, where eps is small and assures that the interplay between t_end, N, dt is correct
• the observer is called after the last step
• the state of the ode x must have the same value as the state in the last call of the observer

[demianmnave]

Here is one reference:

Designing and mapping trimming curves on surfaces using orthogonal projection. Pegna, J and Wolter, F-E. Proc. 16th ASME Design Automation Conf.: Advances in Design Automation, Computer Aided and Computational Design, Vol 1, pp. 235-245, 1990.

If you stick to simple polynomial curves and analytical parametric surfaces (e.g. planes and spheres), the additional coding required to make a test case should be manageable. A somewhat simpler problem to try would be finding the parametric location 't' at which a curve's length is 's'; e.g.:

http://www.geometrictools.com/Documentation/MovingAlongCurveSpecifiedSpeed.pdf

Lines and circles are obvious choices for such a test case since their arc lengths are easy to compute.

Cheers, Demian

[demianmnave]

The main problem is in determining eps, especially since 'state_type' is not restricted to simple floating-point types (as far as I gathered from the documentation). eps would depend upon the accumulated error in the summation, which, at least for IEEE floating-point types, can be computed on-the-fly; see, for example:

Error Estimation of Floating-Point Summation and Dot Product http://www.ti3.tu-harburg.de/paper/rump/Ru11.pdf

Even with this sort of special-purpose error estimation, the chance for surprising results may still remain. Therefore, in my opinion, it is sufficient to document that integrate_const() doesn't integrate at t_end unless dt is an exactly representable multiple of t_end-t (e.g. dt is an IEEE 'double' and an integral multiple of 1/2 for t_end-t=1). Then, having overloads of the 'integrate_n_steps()' function for Controlled Error and Dense Output steppers would satisfy the requirement for precise endpoint integration for all stepper types.

Cheers, Demian

[headmyshoulder]

Yes, this is absolutely reasonable. I want to assure, that the observer is really called in equidistant steps and that the state of the ode after integrate_const is identical to the last state in the observer call, just for consistency.

Btw: At the moment the time is incremented via the += operator. I think should be rethought.

[demianmnave]

That is definitely worthwhile to verify. If you happen to implement a new integerate_n_steps() function, I would be happy to give it a try with my code.

Cheers, Demian

[mariomulansky]

Thanks for the literature hints demian, i will look into this.

I think we should implement integrate_n_steps directly as there we don't have this problem at all and then use this to implement integrate_const. currently the implementation is the other way around. if we do it this way we just have to calculate N in the right way.

[mariomulansky]

i added integrate_n_steps versions for controlled and dense_out steppers, i think we should now use these for implementing the integrate_const steppers. however, i will open a new issue for that.

[demianmnave]

If endpoint integration is important enough to warrant additional implementation, a version of integrate_n_steps() taking (t, t_end, n_steps) can be implemented to exactly include an integration step at t_end.

I definitely agree that integrate_const() could be implemented in terms of the current integrate_n_steps(t, dt, n_steps), although n_steps may have to be computed carefully to avoid an off-by-one error due to roundoff.

[headmyshoulder]

Can we close this issue? I think it is fixed,