headmyshoulder/odeint-v2 — Issue #4: integrate_const with dense output stepper misses final observer call
Author: demianmnave
Created: 2012-02-03 19:58:57+00:00
State: open

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