Merge branch 'master' of github.com:headmyshoulder/odeint-v2 into develop

doc/details_generation_functions.qbk

@@ -18,7 +18,8 @@ In the __tutorial we have learned how we can use the generation functions `make_
 
 [generation_functions_syntax_auto]
 
-The first two parameters are the absolute and the relative error tolerances and the third parameter is the stepper. In C++03 you can infer the type from the `result_of` mechanism:
+The first two parameters are the absolute and the relative error tolerances and the third parameter is the stepper. Additionally, a second version exists where additionally a maximal step size is supplied which ensures the the step size is not increased above this value.
+ In C++03 you can infer the type from the `result_of` mechanism:
 
 [generation_functions_syntax_result_of]
 

doc/details_integrate_functions.qbk

@@ -9,11 +9,15 @@
   http://www.boost.org/LICENSE_1_0.txt)
 =============================================================================/]
 
+[def _max_step_checker_   [classref boost::numeric::odeint::max_step_checker `max_step_checker`]]
+
 [section Integrate functions]
 
 Integrate functions perform the time evolution of a given ODE from some
 starting time ['t[sub 0]] to a given end time ['t[sub 1]] and starting at state ['x[sub 0]] by subsequent calls of a given stepper's `do_step` function.
-Additionally, the user can provide an __observer to analyze the state during time evolution.
+Additionally, the user can provide an __observer to analyze the state during time evolution, and
+a _max_step_checker_ to throw an exception if too many steps are taken between observer calls (i.e. too
+small step size).
 There are five different integrate functions which have different strategies on when to call the observer function during integration.
 All of the integrate functions except `integrate_n_steps` can be called with any stepper following one of the stepper concepts: __stepper , __error_stepper , __controlled_stepper , __dense_output_stepper.
 Depending on the abilities of the stepper, the integrate functions make use of step-size control or dense output.
@@ -28,10 +32,14 @@ We start with explaining `integrate_const`:
 
 `integrate_const( stepper , system , x0 , t0 , t1 , dt , observer )`
 
+`integrate_const( stepper , system , x0 , t0 , t1 , dt , observer , max_step_checker )`
+
 These integrate the ODE given by `system` with subsequent steps from `stepper`.
 Integration start at `t0` and `x0` and ends at some ['t' = t[sub 0] + n dt] with /n/ such that ['t[sub 1] - dt < t' <= t[sub 1]].
 `x0` is changed to the approximative solution ['x(t')] at the end of integration.
 If provided, the `observer` is invoked at times ['t[sub 0]], ['t[sub 0] + dt], ['t[sub 0] + 2dt], ... ,['t'].
+If provided, the `max_step_checker` counts the number of steps between observer calls and throws a
+`no_progress_error` this exceeds some limit (default: 500).
 `integrate_const` returns the number of steps performed during the integration.
 Note that if you are using a simple __stepper or __error_stepper and want to make exactly `n` steps you should prefer the `integrate_n_steps` function below.
 
@@ -54,9 +62,13 @@ steps. The integration is then performed until the time `t0+n*dt`.
 
 `integrate_n_steps( stepper , system , x0 , t0 , dt , n , observer )`
 
+`integrate_n_steps( stepper , system , x0 , t0 , dt , n , observer , max_step_checker )`
+
 Integrates the ODE given by `system` with subsequent steps from `stepper` starting at ['x[sub 0]] and ['t[sub 0]].
 If provided, `observer` is called after every step and at the beginning with
 `t0`, similar as above.
+Again, providing a `max_step_checker` will throw a `no_progress_error` if too many steps are performed
+between observer calls.
 The approximate result for ['x( t[sub 0] + n dt )] is stored in `x0`.
 This function returns the end time `t0 + n*dt`.
 
@@ -76,6 +88,11 @@ Integration start at `t0` and `x0` and ends at ['t[sub 1]].
 If provided, the `observer` is called after each step (and before the first step at `t0`).
 `integrate_adaptive` returns the number of steps performed during the integration.
 
+[note `integrate_adaptive` by design performs an observer call after each time step. Hence
+there is no need for a _max_step_checker_ as only exactly one step is ever performed between
+observer calls.
+]
+
 * If `stepper` is a __stepper or __error_stepper then `dt` is the step size used for integration and `integrate_adaptive` behaves like `integrate_const` except that for the last step the step size is reduced to ensure we end exactly at `t1`.
 If provided, the observer is called at each step.
 * If `stepper` is a __controlled_stepper then `dt` is the initial step size.
@@ -99,6 +116,12 @@ Integration starts at `*times_start` and ends exactly at `*(times_end-1)`.
 `x0` contains the approximate solution at the end point of integration.
 This function requires an observer which is invoked at the subsequent times `*times_start++` until `times_start == times_end`.
 If called with a __boost_range `time_range` the function behaves the same with `times_start = boost::begin( time_range )` and `times_end = boost::end( time_range )`.
+Additionally, a _max_step_checker_ can be provided, e.g.:
+
+`integrate_times( stepper , system , x0 , times_start , times_end , dt , observer , max_step_checker)`
+
+As above, this will throw a `no_progress_error` if too many steps are performed between observer calls.
+
 `integrate_times` returns the number of steps performed during the integration.
 
 * If `stepper` is a __stepper or __error_stepper `dt` is the step size used for integration.

doc/odeint.qbk

@@ -18,7 +18,7 @@
      [id odeint]
      [dirname odeint]
      [authors [Ahnert, Karsten], [Mulansky, Mario]]
-     [copyright 2009-2012 Karsten Ahnert and Mario Mulansky]
+     [copyright 2009-2015 Karsten Ahnert and Mario Mulansky]
      [category math]
      [purpose 
         Numerical integration of ordinary differential equations.

examples/Jamfile.v2

@@ -40,7 +40,8 @@ exe bind_member_functions_cpp11 : bind_member_functions_cpp11.cpp : <cxxflags>-s
 exe molecular_dynamics : molecular_dynamics.cpp : <cxxflags>-std=c++0x ;
 exe molecular_dynamics_cells : molecular_dynamics_cells.cpp : <cxxflags>-std=c++0x ;
 exe abm_precision : abm_precision.cpp ;
-exe integrate_times : integrate_times.cpp ;  
+exe integrate_times : integrate_times.cpp ;
+exe find_crossing : find_crossing.cpp : <cxxflags>-std=c++0x ;
 
 build-project multiprecision ;
 # build-project mtl ;

examples/find_crossing.cpp

@@ -0,0 +1,134 @@
+/*
+ * find_crossing.cpp
+ *
+ * Finds the energy threshold crossing for a damped oscillator.
+ * The algorithm uses a dense out stepper with find_if to first find an
+ * interval containing the threshold crossing and the utilizes the dense out
+ * functionality with a bisection to further refine the interval until some
+ * desired precision is reached.
+ *
+ * Copyright 2015 Mario Mulansky
+ *
+ * Distributed under the Boost Software License, Version 1.0.
+ * (See accompanying file LICENSE_1_0.txt or
+ * copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+
+
+#include <iostream>
+#include <utility>
+#include <algorithm>
+#include <array>
+
+#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
+#include <boost/numeric/odeint/stepper/generation.hpp>
+#include <boost/numeric/odeint/iterator/adaptive_iterator.hpp>
+
+namespace odeint = boost::numeric::odeint;
+
+typedef std::array<double, 2> state_type;
+
+const double gam = 1.0;  // damping strength
+
+void damped_osc(const state_type &x, state_type &dxdt, const double /*t*/)
+{
+    dxdt[0] = x[1];
+    dxdt[1] = -x[0] - gam * x[1];
+}
+
+
+struct energy_condition {
+
+    // defines the threshold crossing in terms of a boolean functor
+
+    double m_min_energy;
+
+    energy_condition(const double min_energy)
+        : m_min_energy(min_energy) { }
+
+    double energy(const state_type &x) {
+        return 0.5 * x[1] * x[1] + 0.5 * x[0] * x[0];
+    }
+
+    bool operator()(const state_type &x) {
+        // becomes true if the energy becomes smaller than the threshold
+        return energy(x) <= m_min_energy;
+    }
+};
+
+
+template<class System, class Condition>
+std::pair<double, state_type>
+find_condition(state_type &x0, System sys, Condition cond,
+               const double t_start, const double t_end, const double dt,
+               const double precision = 1E-6) {
+
+    // integrates an ODE until some threshold is crossed
+    // returns time and state at the point of the threshold crossing
+    // if no threshold crossing is found, some time > t_end is returned
+
+    auto stepper = odeint::make_dense_output(1.0e-6, 1.0e-6,
+                                             odeint::runge_kutta_dopri5<state_type>());
+
+    auto ode_range = odeint::make_adaptive_range(std::ref(stepper), sys, x0,
+                                                 t_start, t_end, dt);
+
+    // find the step where the condition changes
+    auto found_iter = std::find_if(ode_range.first, ode_range.second, cond);
+
+    if(found_iter == ode_range.second)
+    {
+        // no threshold crossing -> return time after t_end and ic
+        return std::make_pair(t_end + dt, x0);
+    }
+
+    // the dense out stepper now covers the interval where the condition changes
+    // improve the solution by bisection
+    double t0 = stepper.previous_time();
+    double t1 = stepper.current_time();
+    double t_m;
+    state_type x_m;
+    // use odeint's resizing functionality to allocate memory for x_m
+    odeint::adjust_size_by_resizeability(x_m, x0,
+                                         typename odeint::is_resizeable<state_type>::type());
+    while(std::abs(t1 - t0) > precision) {
+        t_m = 0.5 * (t0 + t1);  // get the mid point time
+        stepper.calc_state(t_m, x_m); // obtain the corresponding state
+        if (cond(x_m))
+            t1 = t_m;  // condition changer lies before midpoint
+        else
+            t0 = t_m;  // condition changer lies after midpoint
+    }
+    // we found the interval of size eps, take it's midpoint as final guess
+    t_m = 0.5 * (t0 + t1);
+    stepper.calc_state(t_m, x_m);
+    return std::make_pair(t_m, x_m);
+}
+
+
+int main(int argc, char **argv)
+{
+    state_type x0 = {{10.0, 0.0}};
+    const double t_start = 0.0;
+    const double t_end = 10.0;
+    const double dt = 0.1;
+    const double threshold = 0.1;
+
+    energy_condition cond(threshold);
+    state_type x_cond;
+    double t_cond;
+    std::tie(t_cond, x_cond) = find_condition(x0, damped_osc, cond,
+                                              t_start, t_end, dt, 1E-6);
+    if(t_cond > t_end)
+    {
+        // time after t_end -> no threshold crossing within [t_start, t_end]
+        std::cout << "No threshold crossing found." << std::endl;
+    } else
+    {
+        std::cout.precision(16);
+        std::cout << "Time of energy threshold crossing: " << t_cond << std::endl;
+        std::cout << "State: [" << x_cond[0] << " , " << x_cond[1] << "]" << std::endl;
+        std::cout << "Energy: " << cond.energy(x_cond) << std::endl;
+    }
+}

examples/generation_functions.cpp

@@ -64,6 +64,13 @@ struct controller_factory< custom_stepper , custom_controller >
     {
         return custom_controller();
     }
+
+    custom_controller operator()( double abs_tol , double rel_tol , double max_dt ,
+                                  const custom_stepper & ) const
+    {
+        // version with maximal allowed step size max_dt
+        return custom_controller();
+    }
 };
 
 } } }
@@ -77,6 +84,9 @@ int main( int argc , char **argv )
         /*
         //[ generation_functions_syntax_auto
         auto stepper1 = make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() );
+        // or with max step size limit:
+        // auto stepper1 = make_controlled( 1.0e-6 , 1.0e-6 , 0.01, stepper_type() );
+
         auto stepper2 = make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() );
         //]
         */

examples/lorenz_point.cpp

@@ -58,7 +58,7 @@ public:
 // only required for steppers with error control
 point3D operator/( const point3D &p1 , const point3D &p2 )
 {
-    return point3D( p1.x/p2.x , p1.y/p2.y , p1.z/p1.z );
+    return point3D( p1.x/p2.x , p1.y/p2.y , p1.z/p2.z );
 }
 
 point3D abs( const point3D &p )

examples/my_vector.cpp

@@ -14,7 +14,7 @@
 #include <boost/numeric/odeint.hpp>
 
 //[my_vector
-template< int MAX_N >
+template< size_t MAX_N >
 class my_vector
 {
     typedef std::vector< double > vector;
@@ -103,7 +103,10 @@ int main()
     state_type x(3);
     x[0] = 5.0 ; x[1] = 10.0 ; x[2] = 10.0;
 
-    // my_vector works with range_algebra as it implements 
+    // make sure resizing is ON
+    BOOST_STATIC_ASSERT( is_resizeable<state_type>::value == true );
+
+    // my_vector works with range_algebra as it implements
     // the required parts of a container interface
     // no further work is required
 

include/boost/numeric/odeint.hpp

@@ -71,6 +71,7 @@
 #include <boost/numeric/odeint/integrate/integrate_times.hpp>
 
 #include <boost/numeric/odeint/integrate/observer_collection.hpp>
+#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
 
 #include <boost/numeric/odeint/stepper/generation.hpp>
 

include/boost/numeric/odeint/integrate/check_adapter.hpp

@@ -0,0 +1,222 @@
+/*
+ [auto_generated]
+ boost/numeric/odeint/integrate/check_adapter.hpp
+
+ [begin_description]
+ Adapters to add checking facility to stepper and observer
+ [end_description]
+
+ Copyright 2015 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_CHECK_ADAPTER_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_INTEGRATE_CHECK_ADAPTER_HPP_INCLUDED
+
+#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
+#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+template<class Stepper, class Checker,
+         class StepperCategory = typename base_tag<typename Stepper::stepper_category>::type>
+class checked_stepper;
+
+
+/**
+ * \brief Adapter to combine basic stepper and checker.
+ */
+template<class Stepper, class Checker>
+class checked_stepper<Stepper, Checker, stepper_tag>
+{
+
+public:
+    typedef Stepper stepper_type;
+    typedef Checker checker_type;
+    // forward stepper typedefs
+    typedef typename stepper_type::state_type state_type;
+    typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::deriv_type deriv_type;
+    typedef typename stepper_type::time_type time_type;
+
+private:
+    stepper_type &m_stepper;
+    checker_type &m_checker;
+
+public:
+    /**
+     * \brief Construct the checked_stepper.
+     */
+    checked_stepper(stepper_type &stepper, checker_type &checker)
+            : m_stepper(stepper), m_checker(checker) { }
+
+    /**
+     * \brief forward of the do_step method
+     */
+    template<class System, class StateInOut>
+    void do_step(System system, StateInOut &state, const time_type t, const time_type dt)
+    {
+        // do the step
+        m_stepper.do_step(system, state, t, dt);
+        // call the checker
+        m_checker();
+    }
+};
+
+
+/**
+ * \brief Adapter to combine controlled stepper and checker.
+ */
+template<class ControlledStepper, class Checker>
+class checked_stepper<ControlledStepper, Checker, controlled_stepper_tag>
+{
+
+public:
+    typedef ControlledStepper stepper_type;
+    typedef Checker checker_type;
+    // forward stepper typedefs
+    typedef typename stepper_type::state_type state_type;
+    typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::deriv_type deriv_type;
+    typedef typename stepper_type::time_type time_type;
+
+private:
+    stepper_type &m_stepper;
+    checker_type &m_checker;
+
+public:
+    /**
+     * \brief Construct the checked_stepper.
+     */
+    checked_stepper(stepper_type &stepper, checker_type &checker)
+            : m_stepper(stepper), m_checker(checker) { }
+
+    /**
+     * \brief forward of the do_step method
+     */
+    template< class System , class StateInOut >
+    controlled_step_result try_step( System system , StateInOut &state , time_type &t , time_type &dt )
+    {
+        // do the step
+        if( m_stepper.try_step(system, state, t, dt) == success )
+        {
+            // call the checker if step was successful
+            m_checker();
+            return success;
+        } else
+        {
+            // step failed -> return fail
+            return fail;
+        }
+    }
+};
+
+
+/**
+ * \brief Adapter to combine dense out stepper and checker.
+ */
+template<class DenseOutStepper, class Checker>
+class checked_stepper<DenseOutStepper, Checker, dense_output_stepper_tag>
+{
+
+public:
+    typedef DenseOutStepper stepper_type;
+    typedef Checker checker_type;
+    // forward stepper typedefs
+    typedef typename stepper_type::state_type state_type;
+    typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::deriv_type deriv_type;
+    typedef typename stepper_type::time_type time_type;
+
+private:
+    stepper_type &m_stepper;
+    checker_type &m_checker;
+
+public:
+    /**
+     * \brief Construct the checked_stepper.
+     */
+    checked_stepper(stepper_type &stepper, checker_type &checker)
+            : m_stepper(stepper), m_checker(checker) { }
+
+
+    template< class System >
+    std::pair< time_type , time_type > do_step( System system )
+    {
+        m_checker();
+        return m_stepper.do_step(system);
+    }
+
+    /* provide the remaining dense out stepper interface */
+    template< class StateType >
+    void initialize( const StateType &x0 , time_type t0 , time_type dt0 )
+    { m_stepper.initialize(x0, t0, dt0); }
+
+
+    template< class StateOut >
+    void calc_state( time_type t , StateOut &x ) const
+    { m_stepper.calc_state(t, x); }
+
+    template< class StateOut >
+    void calc_state( time_type t , const StateOut &x ) const
+    { m_stepper.calc_state(t, x); }
+
+    const state_type& current_state( void ) const
+    { return m_stepper.current_state(); }
+
+    time_type current_time( void ) const
+    { return m_stepper.current_time(); }
+
+    const state_type& previous_state( void ) const
+    { return m_stepper.previous_state(); }
+
+    time_type previous_time( void ) const
+    { return m_stepper.previous_time(); }
+
+    time_type current_time_step( void ) const
+    { return m_stepper.current_time_step(); }
+
+};
+
+
+/**
+ * \brief Adapter to combine observer and checker.
+ */
+template<class Observer, class Checker>
+class checked_observer
+{
+public:
+    typedef Observer observer_type;
+    typedef Checker checker_type;
+
+private:
+    observer_type &m_observer;
+    checker_type &m_checker;
+
+public:
+    checked_observer(observer_type &observer, checker_type &checker)
+            : m_observer(observer), m_checker(checker)
+    {}
+
+    template< class State , class Time >
+    void operator()(const State& state, Time t) const
+    {
+        // call the observer
+        m_observer(state, t);
+        // reset the checker
+        m_checker.reset();
+    }
+};
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+#endif

include/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp

@@ -7,7 +7,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
 
  Distributed under the Boost Software License, Version 1.0.
@@ -25,6 +25,7 @@
 
 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
+#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
 #include <boost/numeric/odeint/integrate/detail/integrate_const.hpp>
 #include <boost/numeric/odeint/util/bind.hpp>
 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
@@ -41,7 +42,7 @@ namespace odeint {
 namespace detail {
 
 // forward declaration
-template< class Stepper , class System , class State , class Time , class Observer>
+template< class Stepper , class System , class State , class Time , class Observer >
 size_t integrate_const(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time end_time , Time dt ,
@@ -74,7 +75,7 @@ size_t integrate_adaptive(
 
 
 /*
- * classical integrate adaptive
+ * integrate adaptive for controlled stepper
  */
 template< class Stepper , class System , class State , class Time , class Observer >
 size_t integrate_adaptive(
@@ -86,8 +87,7 @@ size_t integrate_adaptive(
     typename odeint::unwrap_reference< Observer >::type &obs = observer;
     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
 
-    const size_t max_attempts = 1000;
-    const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found.";
+    failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
     size_t count = 0;
     while( less_with_sign( start_time , end_time , dt ) )
     {
@@ -97,15 +97,14 @@ size_t integrate_adaptive(
             dt = end_time - start_time;
         }
 
-        size_t trials = 0;
         controlled_step_result res;
         do
         {
             res = st.try_step( system , start_state , start_time , dt );
-            ++trials;
+            fail_checker();  // check number of failed steps
         }
-        while( ( res == fail ) && ( trials < max_attempts ) );
-        if( trials == max_attempts ) BOOST_THROW_EXCEPTION( std::overflow_error( error_string ) );
+        while( res == fail );
+        fail_checker.reset();  // if we reach here, the step was successful -> reset fail checker
 
         ++count;
     }

include/boost/numeric/odeint/integrate/detail/integrate_const.hpp

@@ -6,7 +6,7 @@
  integrate const implementation
  [end_description]
 
- Copyright 2012 Mario Mulansky
+ Copyright 2012-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
  Copyright 2012 Karsten Ahnert
 
@@ -43,12 +43,13 @@ template< class Stepper , class System , class State , class Time , class Observ
 size_t integrate_const(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time end_time , Time dt ,
-        Observer observer , stepper_tag 
+        Observer observer , stepper_tag
 )
 {
+
     typename odeint::unwrap_reference< Observer >::type &obs = observer;
     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
-    
+
     Time time = start_time;
     int step = 0;
     // cast time+dt explicitely in case of expression templates (e.g. multiprecision)
@@ -72,11 +73,11 @@ template< class Stepper , class System , class State , class Time , class Observ
 size_t integrate_const(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time end_time , Time dt ,
-        Observer observer , controlled_stepper_tag 
+        Observer observer , controlled_stepper_tag
 )
 {
     typename odeint::unwrap_reference< Observer >::type &obs = observer;
-    
+
     Time time = start_time;
     const Time time_step = dt;
     int real_steps = 0;
@@ -85,8 +86,10 @@ size_t integrate_const(
     while( less_eq_with_sign( static_cast<Time>(time+time_step) , end_time , dt ) )
     {
         obs( start_state , time );
-        real_steps += detail::integrate_adaptive( stepper , system , start_state , time , time+time_step , dt ,
-                                    null_observer() , controlled_stepper_tag() );
+        // integrate_adaptive_checked uses the given checker to throw if an overflow occurs
+        real_steps += detail::integrate_adaptive(stepper, system, start_state, time,
+                                                 static_cast<Time>(time + time_step), dt,
+                                                 null_observer(), controlled_stepper_tag());
         // direct computation of the time avoids error propagation happening when using time += dt
         // we need clumsy type analysis to get boost units working here
         step++;
@@ -102,12 +105,12 @@ template< class Stepper , class System , class State , class Time , class Observ
 size_t integrate_const(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time end_time , Time dt ,
-        Observer observer , dense_output_stepper_tag 
+        Observer observer , dense_output_stepper_tag
 )
 {
     typename odeint::unwrap_reference< Observer >::type &obs = observer;
     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
-    
+
     Time time = start_time;
     
     st.initialize( start_state , time , dt );
@@ -117,7 +120,7 @@ size_t integrate_const(
     int obs_step( 1 );
     int real_step( 0 );
     
-    while( less_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
+    while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
     {
         while( less_eq_with_sign( time , st.current_time() , dt ) )
         {
@@ -148,6 +151,7 @@ size_t integrate_const(
         
     }
     // last observation, if we are still in observation interval
+    // might happen due to finite precision problems when computing the the time
     if( less_eq_with_sign( time , end_time , dt ) )
     {
         st.calc_state( time , start_state );

include/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp

@@ -6,7 +6,7 @@
  integrate steps implementation
  [end_description]
 
- Copyright 2012 Mario Mulansky
+ Copyright 2012-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
  Copyright 2012 Karsten Ahnert
 
@@ -32,10 +32,10 @@ namespace detail {
 
 // forward declaration
 template< class Stepper , class System , class State , class Time , class Observer >
-size_t integrate_adaptive(
+size_t integrate_adaptive_checked(
         Stepper stepper , System system , State &start_state ,
         Time &start_time , Time end_time , Time &dt ,
-        Observer observer , controlled_stepper_tag
+        Observer observer, controlled_stepper_tag
 );
 
 
@@ -66,7 +66,7 @@ Time integrate_n_steps(
 
 
 /* controlled version */
-template< class Stepper , class System , class State , class Time , class Observer>
+template< class Stepper , class System , class State , class Time , class Observer >
 Time integrate_n_steps(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time dt , size_t num_of_steps ,
@@ -80,8 +80,9 @@ Time integrate_n_steps(
     for( size_t step = 0; step < num_of_steps ; ++step )
     {
         obs( start_state , time );
-        detail::integrate_adaptive( stepper , system , start_state , time , static_cast<Time>(time+time_step) , dt ,
-                null_observer() , controlled_stepper_tag() );
+        // integrate_adaptive_checked uses the given checker to throw if an overflow occurs
+        detail::integrate_adaptive(stepper, system, start_state, time, static_cast<Time>(time + time_step), dt,
+                                   null_observer(), controlled_stepper_tag());
         // direct computation of the time avoids error propagation happening when using time += dt
         // we need clumsy type analysis to get boost units working here
         time = start_time + static_cast< typename unit_value_type<Time>::type >(step+1) * time_step;
@@ -93,7 +94,7 @@ Time integrate_n_steps(
 
 
 /* dense output version */
-template< class Stepper , class System , class State , class Time , class Observer>
+template< class Stepper , class System , class State , class Time , class Observer >
 Time integrate_n_steps(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time dt , size_t num_of_steps ,
@@ -135,6 +136,7 @@ Time integrate_n_steps(
         }
     }
 
+    // make sure we really end exactly where we should end
     while( st.current_time() < end_time )
     {
         if( less_with_sign( end_time ,
@@ -144,7 +146,7 @@ Time integrate_n_steps(
         st.do_step( system );
     }
 
-    // observation at end point, only if we ended exactly on the end-point (or above due to finite precision)
+    // observation at final point
     obs( st.current_state() , end_time );
 
     return time;

include/boost/numeric/odeint/integrate/detail/integrate_times.hpp

@@ -6,7 +6,7 @@
  Default integrate times implementation.
  [end_description]
 
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
  Copyright 2012 Karsten Ahnert
  Copyright 2012 Christoph Koke
 
@@ -26,6 +26,7 @@
 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
+#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
 
 
 namespace boost {
@@ -38,15 +39,18 @@ namespace detail {
 /*
  * integrate_times for simple stepper
  */
-template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
+template<class Stepper, class System, class State, class TimeIterator, class Time, class Observer>
 size_t integrate_times(
         Stepper stepper , System system , State &start_state ,
         TimeIterator start_time , TimeIterator end_time , Time dt ,
         Observer observer , stepper_tag
 )
 {
-    typename odeint::unwrap_reference< Observer >::type &obs = observer;
-    typename odeint::unwrap_reference< Stepper >::type &st = stepper;
+    typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+    typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+
+    stepper_type &st = stepper;
+    observer_type &obs = observer;
     typedef typename unit_value_type<Time>::type time_type;
 
     size_t steps = 0;
@@ -82,12 +86,10 @@ size_t integrate_times(
     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
     typedef typename unit_value_type<Time>::type time_type;
 
-    const size_t max_attempts = 1000;
-    const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found.";
+    failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
     size_t steps = 0;
     while( true )
     {
-        size_t fail_steps = 0;
         Time current_time = *start_time++;
         obs( start_state , current_time );
         if( start_time == end_time )
@@ -99,15 +101,16 @@ size_t integrate_times(
             if( st.try_step( system , start_state , current_time , current_dt ) == success )
             {
                 ++steps;
+                // successful step -> reset the fail counter, see #173
+                fail_checker.reset();
                 // continue with the original step size if dt was reduced due to observation
                 dt = max_abs( dt , current_dt );
             }
             else
             {
-                ++fail_steps;
+                fail_checker();  // check for possible overflow of failed steps in step size adjustment
                 dt = current_dt;
             }
-            if( fail_steps == max_attempts ) BOOST_THROW_EXCEPTION( std::overflow_error( error_string ));
         }
     }
     return steps;
@@ -125,6 +128,7 @@ size_t integrate_times(
 {
     typename odeint::unwrap_reference< Observer >::type &obs = observer;
     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
+
     typedef typename unit_value_type<Time>::type time_type;
 
     if( start_time == end_time )

include/boost/numeric/odeint/integrate/integrate_adaptive.hpp

@@ -7,7 +7,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
 
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or

include/boost/numeric/odeint/integrate/integrate_const.hpp

@@ -8,7 +8,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
 
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or
@@ -23,6 +23,7 @@
 
 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
 #include <boost/numeric/odeint/integrate/null_observer.hpp>
+#include <boost/numeric/odeint/integrate/check_adapter.hpp>
 #include <boost/numeric/odeint/integrate/detail/integrate_const.hpp>
 #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
 
@@ -31,91 +32,120 @@ namespace numeric {
 namespace odeint {
 
 
-
-
-
 /*
  * Integrates with constant time step dt.
  */
-template< class Stepper , class System , class State , class Time , class Observer >
+template<class Stepper, class System, class State, class Time, class Observer, class StepOverflowChecker>
 size_t integrate_const(
-        Stepper stepper , System system , State &start_state ,
-        Time start_time , Time end_time , Time dt ,
-        Observer observer
-)
-{
-    typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
+        Stepper stepper, System system, State &start_state,
+        Time start_time, Time end_time, Time dt,
+        Observer observer, StepOverflowChecker checker
+) {
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
     // we want to get as fast as possible to the end
-    if( boost::is_same< null_observer , Observer >::value )
-    {
+    // no overflow checks needed
+    if (boost::is_same<null_observer, Observer>::value) {
         return detail::integrate_adaptive(
-                stepper , system , start_state ,
-                start_time , end_time  , dt ,
-                observer , stepper_category() );
+                stepper, system, start_state,
+                start_time, end_time, dt,
+                observer, stepper_category());
+    }
+    else {
+        // unwrap references
+        typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+        typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+        typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+
+        return detail::integrate_const(checked_stepper<stepper_type, checker_type>(stepper, checker),
+                                       system, start_state,
+                                       start_time, end_time, dt,
+                                       checked_observer<observer_type, checker_type>(observer, checker),
+                                       stepper_category());
     }
-    else
-    {
-        return detail::integrate_const( stepper , system , start_state , 
-                                        start_time , end_time , dt ,
-                                        observer , stepper_category() );
-      }
 }
 
 /**
- * \brief Second version to solve the forwarding problem, 
- * can be called with Boost.Range as start_state.
- */
-template< class Stepper , class System , class State , class Time , class Observer >
+* \brief Second version to solve the forwarding problem,
+* can be called with Boost.Range as start_state.
+*/
+template<class Stepper, class System, class State, class Time, class Observer, class StepOverflowChecker >
 size_t integrate_const(
-        Stepper stepper , System system , const State &start_state ,
-        Time start_time , Time end_time , Time dt ,
-        Observer observer
-)
-{
-    typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
+        Stepper stepper, System system, const State &start_state,
+        Time start_time, Time end_time, Time dt,
+        Observer observer, StepOverflowChecker checker
+) {
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
     // we want to get as fast as possible to the end
-    if( boost::is_same< null_observer , Observer >::value )
-    {
+
+    if (boost::is_same<null_observer, Observer>::value) {
         return detail::integrate_adaptive(
-                stepper , system , start_state ,
-                start_time , end_time  , dt ,
-                observer , stepper_category() );
+                stepper, system, start_state,
+                start_time, end_time, dt,
+                observer, stepper_category());
     }
-    else
-    {
-        return detail::integrate_const( stepper , system , start_state , 
-                                        start_time , end_time , dt ,
-                                        observer , stepper_category() );
+    else {
+        typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+        typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+        typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+
+        return detail::integrate_const(checked_stepper<stepper_type, checker_type>(stepper, checker),
+                                       system, start_state,
+                                       start_time, end_time, dt,
+                                       checked_observer<observer_type, checker_type>(observer, checker),
+                                       stepper_category());
     }
 }
 
 
-
-
-
 /**
- * \brief integrate_const without observer calls
- */
-template< class Stepper , class System , class State , class Time >
+* \brief integrate_const without step overflow checker
+*/
+template<class Stepper, class System, class State, class Time, class Observer>
 size_t integrate_const(
-        Stepper stepper , System system , State &start_state ,
-        Time start_time , Time end_time , Time dt
-)
+        Stepper stepper, System system, State &start_state,
+        Time start_time, Time end_time, Time dt, Observer observer)
 {
-    return integrate_const( stepper , system , start_state , start_time , end_time , dt , null_observer() );
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
+    return detail::integrate_const(stepper, system, start_state,
+                                   start_time, end_time, dt, observer, stepper_category());
 }
 
 /**
- * \brief Second version to solve the forwarding problem,
- * can be called with Boost.Range as start_state.
- */
-template< class Stepper , class System , class State , class Time >
+* \brief Second version to solve the forwarding problem,
+* can be called with Boost.Range as start_state.
+*/
+template<class Stepper, class System, class State, class Time, class Observer>
 size_t integrate_const(
-        Stepper stepper , System system , const State &start_state ,
-        Time start_time , Time end_time , Time dt
-)
-{
-    return integrate_const( stepper , system , start_state , start_time , end_time , dt , null_observer() );
+        Stepper stepper, System system, const State &start_state,
+        Time start_time, Time end_time, Time dt, Observer observer
+) {
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
+    return detail::integrate_const(stepper, system, start_state,
+                                   start_time, end_time, dt, observer, stepper_category());
+}
+
+
+/**
+* \brief integrate_const without observer calls
+*/
+template<class Stepper, class System, class State, class Time>
+size_t integrate_const(
+        Stepper stepper, System system, State &start_state,
+        Time start_time, Time end_time, Time dt
+) {
+    return integrate_const(stepper, system, start_state, start_time, end_time, dt, null_observer());
+}
+
+/**
+* \brief Second version to solve the forwarding problem,
+* can be called with Boost.Range as start_state.
+*/
+template<class Stepper, class System, class State, class Time>
+size_t integrate_const(
+        Stepper stepper, System system, const State &start_state,
+        Time start_time, Time end_time, Time dt
+) {
+    return integrate_const(stepper, system, start_state, start_time, end_time, dt, null_observer());
 }
 
 
@@ -124,30 +154,37 @@ size_t integrate_const(
 
 
 /********* DOXYGEN *********/
-    /**
-     * \fn integrate_const( Stepper stepper , System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer )
-     * \brief Integrates the ODE with constant step size.
-     *
-     * Integrates the ODE defined by system using the given stepper.
-     * This method ensures that the observer is called at constant intervals dt.
-     * If the Stepper is a normal stepper without step size control, dt is also
-     * used for the numerical scheme. If a ControlledStepper is provided, the 
-     * algorithm might reduce the step size to meet the error bounds, but it is 
-     * ensured that the observer is always called at equidistant time points
-     * t0 + n*dt. If a DenseOutputStepper is used, the step size also may vary
-     * and the dense output is used to call the observer at equidistant time
-     * points.
-     *
-     * \param stepper The stepper to be used for numerical integration.
-     * \param system Function/Functor defining the rhs of the ODE.
-     * \param start_state The initial condition x0.
-     * \param start_time The initial time t0.
-     * \param end_time The final integration time tend.
-     * \param dt The time step between observer calls, _not_ necessarily the 
-     * time step of the integration.
-     * \param observer Function/Functor called at equidistant time intervals.
-     * \return The number of steps performed.
-     */
+/**
+ * \fn integrate_const( Stepper stepper , System system , State &start_state , Time start_time ,
+ *                      Time end_time , Time dt , Observer observer , StepOverflowChecker checker )
+ * \brief Integrates the ODE with constant step size.
+ *
+ * Integrates the ODE defined by system using the given stepper.
+ * This method ensures that the observer is called at constant intervals dt.
+ * If the Stepper is a normal stepper without step size control, dt is also
+ * used for the numerical scheme. If a ControlledStepper is provided, the
+ * algorithm might reduce the step size to meet the error bounds, but it is
+ * ensured that the observer is always called at equidistant time points
+ * t0 + n*dt. If a DenseOutputStepper is used, the step size also may vary
+ * and the dense output is used to call the observer at equidistant time
+ * points.
+ * If a max_step_checker is provided as StepOverflowChecker, a
+ * no_progress_error is thrown if too many steps (default: 500) are performed
+ * without progress, i.e. in between observer calls. If no checker is provided,
+ * no such overflow check is performed.
+ *
+ * \param stepper The stepper to be used for numerical integration.
+ * \param system Function/Functor defining the rhs of the ODE.
+ * \param start_state The initial condition x0.
+ * \param start_time The initial time t0.
+ * \param end_time The final integration time tend.
+ * \param dt The time step between observer calls, _not_ necessarily the
+ * time step of the integration.
+ * \param observer [optional] Function/Functor called at equidistant time intervals.
+ * \param checker [optional] Functor to check for step count overflows, if no
+ * checker is provided, no exception is thrown.
+ * \return The number of steps performed.
+ */
 
 } // namespace odeint
 } // namespace numeric

include/boost/numeric/odeint/integrate/integrate_n_steps.hpp

@@ -7,7 +7,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
 
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or
@@ -34,36 +34,81 @@ namespace odeint {
  *
  * the two overloads are needed in order to solve the forwarding problem
  */
-template< class Stepper , class System , class State , class Time , class Observer>
+template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
 Time integrate_n_steps(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time dt , size_t num_of_steps ,
-        Observer observer )
+        Observer observer , StepOverflowChecker checker )
 {
-    typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
+    // unwrap references
+    typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+    typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+    typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+    typedef typename stepper_type::stepper_category stepper_category;
+
     return detail::integrate_n_steps(
-                stepper , system , start_state ,
+                checked_stepper<stepper_type, checker_type>(stepper, checker),
+                system , start_state ,
                 start_time , dt , num_of_steps ,
-                observer , stepper_category() );
+                checked_observer<observer_type, checker_type>(observer, checker),
+                stepper_category() );
 }
 
 /**
  * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
  */
-template< class Stepper , class System , class State , class Time , class Observer >
+template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
 Time integrate_n_steps(
         Stepper stepper , System system , const State &start_state ,
         Time start_time , Time dt , size_t num_of_steps ,
-        Observer observer )
+        Observer observer , StepOverflowChecker checker )
 {
-    typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
+    typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+    typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+    typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+    typedef typename stepper_type::stepper_category stepper_category;
+
     return detail::integrate_n_steps(
-                 stepper , system , start_state ,
-                 start_time , dt , num_of_steps ,
-                 observer , stepper_category() );
+            checked_stepper<stepper_type, checker_type>(stepper, checker),
+            system , start_state ,
+            start_time , dt , num_of_steps ,
+            checked_observer<observer_type, checker_type>(observer, checker),
+            stepper_category() );
 }
 
 
+/**
+* \brief The same function as above, but without checker.
+*/
+template< class Stepper , class System , class State , class Time , class Observer >
+Time integrate_n_steps(
+        Stepper stepper , System system , State &start_state ,
+        Time start_time , Time dt , size_t num_of_steps , Observer observer )
+{
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
+
+    return detail::integrate_n_steps(
+            stepper , system , start_state ,
+            start_time , dt , num_of_steps ,
+            observer , stepper_category() );
+}
+
+/**
+* \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
+*/
+template< class Stepper , class System , class State , class Time , class Observer >
+Time integrate_n_steps(
+        Stepper stepper , System system , const State &start_state ,
+        Time start_time , Time dt , size_t num_of_steps , Observer observer )
+{
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
+
+    return detail::integrate_n_steps(
+            stepper , system , start_state ,
+            start_time , dt , num_of_steps ,
+            observer , stepper_category() );
+}
+
 /**
  * \brief The same function as above, but without observer calls.
  */
@@ -72,7 +117,8 @@ Time integrate_n_steps(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time dt , size_t num_of_steps )
 {
-    return integrate_n_steps( stepper , system , start_state , start_time , dt , num_of_steps , null_observer() );
+    return integrate_n_steps(stepper, system, start_state, start_time,
+                             dt, num_of_steps, null_observer());
 }
 
 /**
@@ -83,7 +129,8 @@ Time integrate_n_steps(
         Stepper stepper , System system , const State &start_state ,
         Time start_time , Time dt , size_t num_of_steps )
 {
-    return integrate_n_steps( stepper , system , start_state , start_time , dt , num_of_steps , null_observer() );
+    return integrate_n_steps(stepper, system, start_state, start_time,
+                             dt, num_of_steps, null_observer());
 }
 
 
@@ -102,6 +149,11 @@ Time integrate_n_steps(
      * t0 + n*dt. If a DenseOutputStepper is used, the step size also may vary
      * and the dense output is used to call the observer at equidistant time
      * points. The final integration time is always t0 + num_of_steps*dt.
+     * If a max_step_checker is provided as StepOverflowChecker, a
+     * no_progress_errror is thrown if too many steps (default: 500) are
+     * performed without progress, i.e. in between observer calls. If no
+     * checker is provided, no such overflow check is performed.
+
      *
      * \param stepper The stepper to be used for numerical integration.
      * \param system Function/Functor defining the rhs of the ODE.
@@ -111,6 +163,8 @@ Time integrate_n_steps(
      * time step of the integration.
      * \param num_of_steps Number of steps to be performed
      * \param observer Function/Functor called at equidistant time intervals.
+     * \param checker [optional] Functor to check for step count overflows, if no
+     * checker is provided, no exception is thrown.
      * \return The number of steps performed.
      */
 

include/boost/numeric/odeint/integrate/integrate_times.hpp

@@ -7,7 +7,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
 
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or
@@ -24,6 +24,7 @@
 
 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
 #include <boost/numeric/odeint/integrate/null_observer.hpp>
+#include <boost/numeric/odeint/integrate/check_adapter.hpp>
 #include <boost/numeric/odeint/integrate/detail/integrate_times.hpp>
 
 namespace boost {
@@ -32,8 +33,91 @@ namespace odeint {
 
 
 /*
+ * \brief Integrates while calling the observer at the time points given by sequence [times_start, time_end)
  * the two overloads are needed in order to solve the forwarding problem
  */
+template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer , class StepOverflowChecker >
+size_t integrate_times(
+        Stepper stepper , System system , State &start_state ,
+        TimeIterator times_start , TimeIterator times_end , Time dt ,
+        Observer observer , StepOverflowChecker checker )
+{
+    // unwrap references
+    typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+    typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+    typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+    typedef typename stepper_type::stepper_category stepper_category;
+
+    // pass on checked stepper and observer
+    // checked_stepper/observer use references internally, so passing by value is fine
+    return detail::integrate_times(
+            checked_stepper<stepper_type, checker_type>(stepper, checker) ,
+            system , start_state ,
+            times_start , times_end , dt ,
+            checked_observer<observer_type, checker_type>(observer, checker),
+            stepper_category() );
+}
+
+/**
+ * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
+ */
+template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer , class StepOverflowChecker >
+size_t integrate_times(
+        Stepper stepper , System system , const State &start_state ,
+        TimeIterator times_start , TimeIterator times_end , Time dt ,
+        Observer observer , StepOverflowChecker checker )
+{
+    typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+    typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+    typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+    typedef typename stepper_type::stepper_category stepper_category;
+
+    stepper_type &st = stepper;
+    observer_type &obs = observer;
+    checker_type &chk = checker;
+
+    return detail::integrate_times(
+            checked_stepper<stepper_type, checker_type>(stepper, checker) ,
+            system , start_state ,
+            times_start , times_end , dt ,
+            checked_observer<observer_type, checker_type>(observer, checker),
+            stepper_category() );
+}
+
+/**
+ * \brief The same function as above, but with the observation times given as range.
+ */
+template< class Stepper , class System , class State , class TimeRange , class Time , class Observer , class StepOverflowChecker >
+size_t integrate_times(
+        Stepper stepper , System system , State &start_state ,
+        const TimeRange &times , Time dt ,
+        Observer observer , StepOverflowChecker checker )
+{
+    return integrate_times(
+            stepper , system , start_state ,
+            boost::begin( times ) , boost::end( times ) , dt , observer , checker );
+}
+
+/**
+ * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
+ */
+template< class Stepper , class System , class State , class TimeRange , class Time , class Observer , class StepOverflowChecker >
+size_t integrate_times(
+        Stepper stepper , System system , const State &start_state ,
+        const TimeRange &times , Time dt ,
+        Observer observer , StepOverflowChecker checker )
+{
+    return integrate_times(
+            stepper , system , start_state ,
+            boost::begin( times ) , boost::end( times ) , dt , observer , checker );
+}
+
+
+
+
+/*
+* The same functions as above, but without a StepOverflowChecker
+*/
 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
 size_t integrate_times(
         Stepper stepper , System system , State &start_state ,
@@ -41,6 +125,7 @@ size_t integrate_times(
         Observer observer )
 {
     typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
+    // simply don't use checked_* adapters
     return detail::integrate_times(
             stepper , system , start_state ,
             times_start , times_end , dt ,
@@ -48,8 +133,8 @@ size_t integrate_times(
 }
 
 /**
- * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
- */
+* \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
+*/
 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
 size_t integrate_times(
         Stepper stepper , System system , const State &start_state ,
@@ -64,8 +149,8 @@ size_t integrate_times(
 }
 
 /**
- * \brief The same function as above, but without observer calls.
- */
+* \brief The same function as above, but with the observation times given as range.
+*/
 template< class Stepper , class System , class State , class TimeRange , class Time , class Observer >
 size_t integrate_times(
         Stepper stepper , System system , State &start_state ,
@@ -78,8 +163,8 @@ size_t integrate_times(
 }
 
 /**
- * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
- */
+* \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
+*/
 template< class Stepper , class System , class State , class TimeRange , class Time , class Observer >
 size_t integrate_times(
         Stepper stepper , System system , const State &start_state ,
@@ -88,12 +173,10 @@ size_t integrate_times(
 {
     return integrate_times(
             stepper , system , start_state ,
-            boost::begin( times ) , boost::end( times ) , dt , observer );
+            boost::begin( times ) , boost::end( times ) , dt , observer);
 }
 
 
-
-
 /********* DOXYGEN ***********/
 
     /**
@@ -110,6 +193,10 @@ size_t integrate_times(
      * If a DenseOutputStepper is provided, the dense output functionality is
      * used to call the observer at the given times. The end time of the 
      * integration is always *(end_time-1).
+     * If a max_step_checker is provided as StepOverflowChecker, a
+     * no_progress_error is thrown if too many steps (default: 500) are
+     * performed without progress, i.e. in between observer calls. If no
+     * checker is provided, no such overflow check is performed.
      *
      * \param stepper The stepper to be used for numerical integration.
      * \param system Function/Functor defining the rhs of the ODE.
@@ -119,6 +206,8 @@ size_t integrate_times(
      * \param dt The time step between observer calls, _not_ necessarily the 
      * time step of the integration.
      * \param observer Function/Functor called at equidistant time intervals.
+     * \param checker [optional] Functor to check for step count overflows, if no
+     * checker is provided, no exception is thrown.
      * \return The number of steps performed.
      */
 

include/boost/numeric/odeint/integrate/max_step_checker.hpp

@@ -0,0 +1,114 @@
+/*
+ [auto_generated]
+ boost/numeric/odeint/integrate/max_step_checker.hpp
+
+ [begin_description]
+ Throws exception if too many steps are performed.
+ [end_description]
+
+ Copyright 2015 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+
+#ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_MAX_STEP_CHECKER_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_INTEGRATE_MAX_STEP_CHECKER_HPP_INCLUDED
+
+#include <stdexcept>
+#include <cstdio>
+
+#include <boost/throw_exception.hpp>
+#include <boost/numeric/odeint/util/odeint_error.hpp>
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+/**
+ * \brief A class for performing overflow checks on the step count in integrate functions.
+ *
+ * Provide an instance of this class to integrate functions if you want to throw a runtime error if
+ * too many steps are performed without progress during the integrate routine.
+ */
+class max_step_checker
+{
+public:
+
+protected:
+    const int m_max_steps;
+    int m_steps;
+
+public:
+    /**
+     * \brief Construct the max_step_checker.
+     * max_steps is the maximal number of iterations allowed before runtime_error is thrown.
+     */
+    max_step_checker(const int max_steps = 500)
+        : m_max_steps(max_steps)
+    {
+        reset();
+    }
+
+    /**
+     * \brief Resets the max_step_checker by setting the internal counter to 0.
+     */
+    void reset()
+    {
+        m_steps = 0;
+    }
+
+    /**
+     * \brief Increases the counter and performs the iteration check
+     */
+    void operator()(void)
+    {
+        if( m_steps++ >= m_max_steps )
+        {
+            char error_msg[200];
+            sprintf(error_msg, "Max number of iterations exceeded (%d).", m_max_steps);
+            BOOST_THROW_EXCEPTION( no_progress_error(error_msg) );
+        }
+    }
+};
+
+
+/**
+ * \brief A class for performing overflow checks on the failed step count in step size adjustments.
+ *
+ * Used internally within the dense output stepper and integrate routines.
+ */
+class failed_step_checker : public max_step_checker
+{
+
+public:
+    /**
+     * \brief Construct the failed_step_checker.
+     * max_steps is the maximal number of iterations allowed before runtime_error is thrown.
+     */
+    failed_step_checker(const int max_steps = 500)
+            : max_step_checker(max_steps)
+    {}
+
+    /**
+     * \brief Increases the counter and performs the iteration check
+     */
+    void operator()(void)
+    {
+        if( m_steps++ >= m_max_steps )
+        {
+            char error_msg[200];
+            sprintf(error_msg, "Max number of iterations exceeded (%d). A new step size was not found.", m_max_steps);
+            BOOST_THROW_EXCEPTION( step_adjustment_error(error_msg) );
+        }
+    }
+};
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+#endif

include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp

@@ -90,9 +90,11 @@ public:
 
     bulirsch_stoer(
         value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
-        value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 )
+        value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 ,
+        time_type max_dt = static_cast<time_type>(0))
         : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , m_midpoint() ,
           m_last_step_rejected( false ) , m_first( true ) ,
+          m_max_dt(max_dt) ,
           m_interval_sequence( m_k_max+1 ) ,
           m_coeff( m_k_max+1 ) ,
           m_cost( m_k_max+1 ) ,
@@ -189,6 +191,14 @@ public:
     template< class System , class StateIn , class DerivIn , class StateOut >
     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
     {
+        if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
+        {
+            // given step size is bigger then max_dt
+            // set limit and return fail
+            dt = m_max_dt;
+            return fail;
+        }
+
         BOOST_USING_STD_MIN();
         BOOST_USING_STD_MAX();
 
@@ -311,6 +321,11 @@ public:
 
         if( !m_last_step_rejected || boost::numeric::odeint::detail::less_with_sign(new_h, dt, dt) )
         {
+            // limit step size
+            if( m_max_dt != static_cast<time_type>(0) )
+            {
+                new_h = detail::min_abs(m_max_dt, new_h);
+            }
             m_dt_last = new_h;
             dt = new_h;
         }
@@ -474,6 +489,7 @@ private:
 
     time_type m_dt_last;
     time_type m_t_last;
+    time_type m_max_dt;
 
     size_t m_current_k_opt;
 
@@ -493,7 +509,7 @@ private:
 
     state_table_type m_table; // sequence of states for extrapolation
 
-    const value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
+    value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
 };
 
 

include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp

@@ -6,7 +6,7 @@
  Implementaiton of the Burlish-Stoer method with dense output
  [end_description]
 
- Copyright 2011-2013 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
  Copyright 2011-2013 Karsten Ahnert
  Copyright 2012 Christoph Koke
 
@@ -43,6 +43,8 @@
 #include <boost/numeric/odeint/util/resizer.hpp>
 #include <boost/numeric/odeint/util/unit_helper.hpp>
 
+#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
+
 #include <boost/type_traits.hpp>
 
 
@@ -97,8 +99,10 @@ public:
     bulirsch_stoer_dense_out(
         value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
         value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 ,
+        time_type max_dt = static_cast<time_type>(0) ,
         bool control_interpolation = false )
-        : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , 
+        : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) ,
+          m_max_dt(max_dt) ,
           m_control_interpolation( control_interpolation) ,
           m_last_step_rejected( false ) , m_first( true ) ,
           m_current_state_x1( true ) ,
@@ -149,6 +153,14 @@ public:
     template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , DerivOut &dxdt_new , time_type &dt )
     {
+        if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
+        {
+            // given step size is bigger then max_dt
+            // set limit and return fail
+            dt = m_max_dt;
+            return fail;
+        }
+
         BOOST_USING_STD_MIN();
         BOOST_USING_STD_MAX();
         using std::pow;
@@ -275,7 +287,14 @@ public:
         }
         //set next stepsize
         if( !m_last_step_rejected || (new_h < dt) )
+        {
+            // limit step size
+            if( m_max_dt != static_cast<time_type>(0) )
+            {
+                new_h = detail::min_abs(m_max_dt, new_h);
+            }
             dt = new_h;
+        }
 
         m_last_step_rejected = reject;
         if( reject )
@@ -301,23 +320,20 @@ public:
     template< class System >
     std::pair< time_type , time_type > do_step( System system )
     {
-        const size_t max_count = 1000;
-
         if( m_first )
         {
             typename odeint::unwrap_reference< System >::type &sys = system;
             sys( get_current_state() , get_current_deriv() , m_t );
         }
 
+        failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
         controlled_step_result res = fail;
         m_t_last = m_t;
-        size_t count = 0;
         while( res == fail )
         {
             res = try_step( system , get_current_state() , get_current_deriv() , m_t , get_old_state() , get_old_deriv() , m_dt );
             m_first = false;
-            if( count++ == max_count )
-                throw std::overflow_error( "bulirsch_stoer : too much iterations!");
+            fail_checker();  // check for overflow of failed steps
         }
         toggle_current_state();
         return std::make_pair( m_t_last , m_t );
@@ -412,15 +428,15 @@ private:
         BOOST_USING_STD_MAX();
         using std::pow;
 
-        value_type expo=1.0/(m_interval_sequence[k-1]);
+        value_type expo = static_cast<value_type>(1)/(m_interval_sequence[k-1]);
         value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
         value_type fac;
         if (error == 0.0)
-            fac=1.0/facmin;
+            fac = static_cast<value_type>(1)/facmin;
         else
         {
             fac = STEPFAC2 / pow BOOST_PREVENT_MACRO_SUBSTITUTION( error / STEPFAC1 , expo );
-            fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( facmin/STEPFAC4 , min BOOST_PREVENT_MACRO_SUBSTITUTION( 1.0/facmin , fac ) );
+            fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( facmin/STEPFAC4 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>(static_cast<value_type>(1)/facmin) , fac ) );
         }
         return h*fac;
     }
@@ -646,6 +662,8 @@ private:
     default_error_checker< value_type, algebra_type , operations_type > m_error_checker;
     modified_midpoint_dense_out< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint;
 
+    time_type m_max_dt;
+
     bool m_control_interpolation;
 
     bool m_last_step_rejected;
@@ -684,7 +702,7 @@ private:
 
     //wrapped_state_type m_a1 , m_a2 , m_a3 , m_a4;
 
-    const value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
+    value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
 };
 
 

include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp

@@ -6,7 +6,7 @@
  [end_description]
 
  Copyright 2010-2013 Karsten Ahnert
- Copyright 2010-2013 Mario Mulansky
+ Copyright 2010-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
 
  Distributed under the Boost Software License, Version 1.0.
@@ -33,6 +33,7 @@
 #include <boost/numeric/odeint/util/state_wrapper.hpp>
 #include <boost/numeric/odeint/util/is_resizeable.hpp>
 #include <boost/numeric/odeint/util/resizer.hpp>
+#include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
 
 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
 #include <boost/numeric/odeint/algebra/default_operations.hpp>
@@ -64,23 +65,24 @@ public:
             value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
             value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
             value_type a_x = static_cast< value_type >( 1 ) ,
-            value_type a_dxdt = static_cast< value_type >( 1 ) )
-    : m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
+            value_type a_dxdt = static_cast< value_type >( 1 ))
+        : m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
     { }
 
 
-    template< class State , class Deriv , class Err , class Time >
+    template< class State , class Deriv , class Err, class Time >
     value_type error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
     {
         return error( algebra_type() , x_old , dxdt_old , x_err , dt );
     }
 
-    template< class State , class Deriv , class Err , class Time >
+    template< class State , class Deriv , class Err, class Time >
     value_type error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
     {
+        using std::abs;
         // this overwrites x_err !
         algebra.for_each3( x_err , x_old , dxdt_old ,
-                typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * get_unit_value( dt ) ) );
+                typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * abs(get_unit_value( dt )) ) );
 
         // value_type res = algebra.reduce( x_err ,
         //        typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0 ) );
@@ -97,9 +99,73 @@ private:
 };
 
 
+template< typename Value, typename Time >
+class default_step_adjuster
+{
+public:
+    typedef Time time_type;
+    typedef Value value_type;
+
+    default_step_adjuster(const time_type max_dt=static_cast<time_type>(0))
+            : m_max_dt(max_dt)
+    {}
 
 
+    time_type decrease_step(time_type dt, const value_type error, const int error_order) const
+    {
+        // returns the decreased time step
+        BOOST_USING_STD_MIN();
+        BOOST_USING_STD_MAX();
+        using std::pow;
 
+        dt *= max
+        BOOST_PREVENT_MACRO_SUBSTITUTION(
+                static_cast<value_type>( static_cast<value_type>(9) / static_cast<value_type>(10) *
+                                         pow(error, static_cast<value_type>(-1) / (error_order - 1))),
+                static_cast<value_type>( static_cast<value_type>(1) / static_cast<value_type> (5)));
+        if(m_max_dt != static_cast<time_type >(0))
+            // limit to maximal stepsize even when decreasing
+            dt = detail::min_abs(dt, m_max_dt);
+        return dt;
+    }
+
+    time_type increase_step(time_type dt, value_type error, const int stepper_order) const
+    {
+        // returns the increased time step
+        BOOST_USING_STD_MIN();
+        BOOST_USING_STD_MAX();
+        using std::pow;
+
+        // adjust the size if dt is smaller than max_dt (providede max_dt is not zero)
+        if(error < 0.5)
+        {
+            // error should be > 0
+            error = max BOOST_PREVENT_MACRO_SUBSTITUTION (
+                    static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(stepper_order) ) ) ,
+                    error);
+            time_type dt_old = dt;
+            //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
+            dt *= static_cast<value_type>(9)/static_cast<value_type>(10) *
+                  pow(error, static_cast<value_type>(-1) / stepper_order);
+            if(m_max_dt != static_cast<time_type >(0))
+                // limit to maximal stepsize
+                dt = detail::min_abs(dt, m_max_dt);
+        }
+        return dt;
+    }
+
+    bool check_step_size_limit(const time_type dt)
+    {
+        if(m_max_dt != static_cast<time_type >(0))
+            return detail::less_eq_with_sign(dt, m_max_dt, dt);
+        return true;
+    }
+
+    time_type get_max_dt() { return m_max_dt; }
+
+private:
+    time_type m_max_dt;
+};
 
 
 
@@ -109,8 +175,10 @@ private:
 template<
 class ErrorStepper ,
 class ErrorChecker = default_error_checker< typename ErrorStepper::value_type ,
-typename ErrorStepper::algebra_type ,
-typename ErrorStepper::operations_type > ,
+    typename ErrorStepper::algebra_type ,
+    typename ErrorStepper::operations_type > ,
+class StepAdjuster = default_step_adjuster< typename ErrorStepper::value_type ,
+    typename ErrorStepper::time_type > ,
 class Resizer = typename ErrorStepper::resizer_type ,
 class ErrorStepperCategory = typename ErrorStepper::stepper_category
 >
@@ -139,11 +207,13 @@ class controlled_runge_kutta ;
  * \tparam Resizer The resizer policy type.
  */
 template<
-class ErrorStepper ,
-class ErrorChecker ,
+class ErrorStepper,
+class ErrorChecker,
+class StepAdjuster,
 class Resizer
 >
-class controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag >
+class controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster, Resizer ,
+        explicit_error_stepper_tag >
 {
 
 public:
@@ -157,13 +227,15 @@ public:
     typedef typename stepper_type::operations_type operations_type;
     typedef Resizer resizer_type;
     typedef ErrorChecker error_checker_type;
+    typedef StepAdjuster step_adjuster_type;
     typedef explicit_controlled_stepper_tag stepper_category;
 
 #ifndef DOXYGEN_SKIP
     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
 
-    typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
+    typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster ,
+            Resizer , explicit_error_stepper_tag > controlled_stepper_type;
 #endif //DOXYGEN_SKIP
 
 
@@ -174,9 +246,10 @@ public:
      */
     controlled_runge_kutta(
             const error_checker_type &error_checker = error_checker_type( ) ,
+            const step_adjuster_type &step_adjuster = step_adjuster_type() ,
             const stepper_type &stepper = stepper_type( )
     )
-    : m_stepper( stepper ) , m_error_checker( error_checker )
+        : m_stepper(stepper), m_error_checker(error_checker) , m_step_adjuster(step_adjuster)
     { }
 
 
@@ -338,58 +411,34 @@ public:
     template< class System , class StateIn , class DerivIn , class StateOut >
     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
     {
-        BOOST_USING_STD_MIN();
-        BOOST_USING_STD_MAX();
-        using std::pow;
+        if( !m_step_adjuster.check_step_size_limit(dt) )
+        {
+            // given dt was above step size limit - adjust and return fail;
+            dt = m_step_adjuster.get_max_dt();
+            return fail;
+        }
 
         m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
 
         // do one step with error calculation
         m_stepper.do_step( system , in , dxdt , t , out , dt , m_xerr.m_v );
 
-        m_max_rel_error = m_error_checker.error( m_stepper.algebra() , in , dxdt , m_xerr.m_v , dt );
+        value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt , m_xerr.m_v , dt );
 
-        if( m_max_rel_error > 1.0 )
+        if( max_rel_err > 1.0 )
         {
-            // error too large - decrease dt ,limit scaling factor to 0.2 and reset state
-            dt *= max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) *
-                                                         pow( m_max_rel_error , static_cast<value_type>(-1) / ( m_stepper.error_order() - 1 ) ) ) ,
-                                                         static_cast<value_type>( static_cast<value_type>(1)/static_cast<value_type> (5) ) );
+            // error too big, decrease step size and reject this step
+            dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
             return fail;
-        }
-        else
+        } else
         {
-            if( m_max_rel_error < 0.5 )
-            {
-                // error should be > 0
-                m_max_rel_error = max BOOST_PREVENT_MACRO_SUBSTITUTION (
-                            static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(m_stepper.stepper_order()) ) ) ,
-                            m_max_rel_error );
-                //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
-                t += dt;
-                dt *= static_cast<value_type>(9)/static_cast<value_type>(10) * pow( m_max_rel_error ,
-                                                                                    static_cast<value_type>(-1) / m_stepper.stepper_order() );
-                return success;
-            }
-            else
-            {
-                t += dt;
-                return success;
-            }
+            // otherwise, increase step size and accept
+            t += dt;
+            dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
+            return success;
         }
     }
 
-    /**
-     * \brief Returns the error of the last step.
-     *
-     * returns The last error of the step.
-     */
-    value_type last_error( void ) const
-    {
-        return m_max_rel_error;
-    }
-
-
     /**
      * \brief Adjust the size of all temporaries in the stepper manually.
      * \param x A state from which the size of the temporaries to be resized is deduced.
@@ -455,6 +504,7 @@ private:
 
     stepper_type m_stepper;
     error_checker_type m_error_checker;
+    step_adjuster_type m_step_adjuster;
 
     resizer_type m_dxdt_resizer;
     resizer_type m_xerr_resizer;
@@ -463,7 +513,6 @@ private:
     wrapped_deriv_type m_dxdt;
     wrapped_state_type m_xerr;
     wrapped_state_type m_xnew;
-    value_type m_max_rel_error;
 };
 
 
@@ -498,9 +547,10 @@ private:
 template<
 class ErrorStepper ,
 class ErrorChecker ,
+class StepAdjuster ,
 class Resizer
 >
-class controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_fsal_tag >
+class controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster , Resizer , explicit_error_stepper_fsal_tag >
 {
 
 public:
@@ -514,13 +564,14 @@ public:
     typedef typename stepper_type::operations_type operations_type;
     typedef Resizer resizer_type;
     typedef ErrorChecker error_checker_type;
+    typedef StepAdjuster step_adjuster_type;
     typedef explicit_controlled_stepper_fsal_tag stepper_category;
 
 #ifndef DOXYGEN_SKIP
     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
 
-    typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
+    typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
 #endif // DOXYGEN_SKIP
 
     /**
@@ -530,9 +581,10 @@ public:
      */
     controlled_runge_kutta(
             const error_checker_type &error_checker = error_checker_type() ,
+            const step_adjuster_type &step_adjuster = step_adjuster_type() ,
             const stepper_type &stepper = stepper_type()
     )
-    : m_stepper( stepper ) , m_error_checker( error_checker ) ,
+    : m_stepper( stepper ) , m_error_checker( error_checker ) , m_step_adjuster(step_adjuster) ,
       m_first_call( true )
     { }
 
@@ -699,10 +751,12 @@ public:
     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t ,
             StateOut &out , DerivOut &dxdt_out , time_type &dt )
     {
-        BOOST_USING_STD_MIN();
-        BOOST_USING_STD_MAX();
-
-        using std::pow;
+        if( !m_step_adjuster.check_step_size_limit(dt) )
+        {
+            // given dt was above step size limit - adjust and return fail;
+            dt = m_step_adjuster.get_max_dt();
+            return fail;
+        }
 
         m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
 
@@ -715,29 +769,14 @@ public:
 
         if( max_rel_err > 1.0 )
         {
-            // error too large - decrease dt ,limit scaling factor to 0.2 and reset state
-            dt *= max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) *
-                                                                 pow( max_rel_err , static_cast<value_type>(-1) / ( m_stepper.error_order() - 1 ) ) ) ,
-                                                         static_cast<value_type>( static_cast<value_type>(1)/static_cast<value_type> (5)) );
+            // error too big, decrease step size and reject this step
+            dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
             return fail;
         }
-        else
-        {
-            if( max_rel_err < 0.5 )
-            {                //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
-                // error should be > 0
-                max_rel_err = max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(m_stepper.stepper_order()) ) ) ,
-                                                                     max_rel_err );
-                t += dt;
-                dt *= static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) * pow( max_rel_err , static_cast<value_type>(-1) / m_stepper.stepper_order() ) );
-                return success;
-            }
-            else
-            {
-                t += dt;
-                return success;
-            }
-        }
+        // otherwise, increase step size and accept
+        t += dt;
+        dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
+        return success;
     }
 
 
@@ -863,6 +902,7 @@ private:
 
     stepper_type m_stepper;
     error_checker_type m_error_checker;
+    step_adjuster_type m_step_adjuster;
 
     resizer_type m_dxdt_resizer;
     resizer_type m_xerr_resizer;
@@ -890,12 +930,14 @@ private:
  * It is used by the controlled_runge_kutta steppers.
  *
  * \tparam Value The value type.
+ * \tparam Time The time type.
  * \tparam Algebra The algebra type.
  * \tparam Operations The operations type.
  */
 
     /**
-     * \fn default_error_checker( value_type eps_abs , value_type eps_rel , value_type a_x , value_type a_dxdt )
+     * \fn default_error_checker( value_type eps_abs , value_type eps_rel , value_type a_x , value_type a_dxdt ,
+     * time_type max_dt)
      * \brief Constructs the error checker.
      *
      * The error is calculated as follows: ???? 
@@ -904,10 +946,11 @@ private:
      * \param eps_rel Relative tolerance level.
      * \param a_x Factor for the weight of the state.
      * \param a_dxdt Factor for the weight of the derivative.
+     * \param max_dt Maximum allowed step size.
      */
     
     /**
-     * \fn error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
+     * \fn error( const State &x_old , const Deriv &dxdt_old , Err &x_err , time_type dt ) const
      * \brief Calculates the error level.
      *
      * If the returned error level is greater than 1, the estimated error was
@@ -922,7 +965,7 @@ private:
      */
 
     /**
-     * \fn error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
+     * \fn error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , time_type dt ) const
      * \brief Calculates the error level using a given algebra.
      *
      * If the returned error level is greater than 1, the estimated error was
@@ -937,6 +980,31 @@ private:
      * \return error
      */
 
+    /**
+     * \fn time_type decrease_step(const time_type dt, const value_type error, const int error_order)
+     * \brief Returns a decreased step size based on the given error and order
+     *
+     * Calculates a new smaller step size based on the given error and its order.
+     *
+     * \param dt The old step size.
+     * \param error The computed error estimate.
+     * \param error_order The error order of the stepper.
+     * \return dt_new The new, reduced step size.
+     */
+
+    /**
+     * \fn time_type increase_step(const time_type dt, const value_type error, const int error_order)
+     * \brief Returns an increased step size based on the given error and order.
+     *
+     * Calculates a new bigger step size based on the given error and its order. If max_dt != 0, the
+     * new step size is limited to max_dt.
+     *
+     * \param dt The old step size.
+     * \param error The computed error estimate.
+     * \param error_order The order of the stepper.
+     * \return dt_new The new, increased step size.
+     */
+
 
 } // odeint
 } // numeric

include/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp

@@ -8,7 +8,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
 
  Distributed under the Boost Software License, Version 1.0.
@@ -37,6 +37,8 @@
 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
 
+#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
+
 namespace boost {
 namespace numeric {
 namespace odeint {
@@ -206,6 +208,15 @@ public:
         return m_t_old;
     }
 
+    /**
+     * \brief Returns the current time step.
+     * \return dt.
+     */
+    time_type current_time_step( void ) const
+    {
+        return m_dt;
+    }
+
 
 private:
 
@@ -312,8 +323,6 @@ public:
     template< class System >
     std::pair< time_type , time_type > do_step( System system )
     {
-        const size_t max_count = 1000;
-
         if( !m_is_deriv_initialized )
         {
             typename odeint::unwrap_reference< System >::type &sys = system;
@@ -321,15 +330,14 @@ public:
             m_is_deriv_initialized = true;
         }
 
+        failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
         controlled_step_result res = fail;
         m_t_old = m_t;
-        size_t count = 0;
         do
         {
             res = m_stepper.try_step( system , get_current_state() , get_current_deriv() , m_t ,
                                       get_old_state() , get_old_deriv() , m_dt );
-            if( count++ == max_count )
-                BOOST_THROW_EXCEPTION( std::overflow_error( "dense_output_controlled_explicit : too much iterations!") );
+            fail_checker();  // check for overflow of failed steps
         }
         while( res == fail );
         toggle_current_state();

include/boost/numeric/odeint/stepper/generation/generation_controlled_runge_kutta.hpp

@@ -34,17 +34,25 @@ struct controller_factory< Stepper , controlled_runge_kutta< Stepper > >
     typedef Stepper stepper_type;
     typedef controlled_runge_kutta< stepper_type > controller_type;
     typedef typename controller_type::error_checker_type error_checker_type;
+    typedef typename controller_type::step_adjuster_type step_adjuster_type;
     typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::value_type time_type;
 
     controller_type operator()( value_type abs_error , value_type rel_error , const stepper_type &stepper )
     {
-        return controller_type( error_checker_type( abs_error , rel_error ) , stepper );
+        return controller_type( error_checker_type( abs_error , rel_error ) ,
+                                step_adjuster_type() , stepper );
+    }
+
+    controller_type operator()( value_type abs_error , value_type rel_error ,
+                                time_type max_dt, const stepper_type &stepper )
+    {
+        return controller_type( error_checker_type( abs_error , rel_error ) ,
+                                step_adjuster_type(max_dt) , stepper );
     }
 };
 
 
-
-
 } // odeint
 } // numeric
 } // boost

include/boost/numeric/odeint/stepper/generation/generation_dense_output_runge_kutta.hpp

@@ -33,12 +33,23 @@ struct dense_output_factory< Stepper , dense_output_runge_kutta< controlled_rung
     typedef Stepper stepper_type;
     typedef controlled_runge_kutta< stepper_type > controller_type;
     typedef typename controller_type::error_checker_type error_checker_type;
+    typedef typename controller_type::step_adjuster_type step_adjuster_type;
     typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::time_type time_type;
     typedef dense_output_runge_kutta< controller_type > dense_output_type;
 
     dense_output_type operator()( value_type abs_error , value_type rel_error , const stepper_type &stepper )
     {
-        return dense_output_type( controller_type( error_checker_type( abs_error , rel_error ) , stepper ) );
+        return dense_output_type( controller_type( error_checker_type( abs_error , rel_error ) ,
+                                                   step_adjuster_type() , stepper ) );
+    }
+
+    dense_output_type operator()( value_type abs_error , value_type rel_error ,
+                                  time_type max_dt , const stepper_type &stepper )
+    {
+        return dense_output_type(
+                controller_type( error_checker_type( abs_error , rel_error) ,
+                                 step_adjuster_type( max_dt ) , stepper ) );
     }
 };
 

include/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp

@@ -54,12 +54,19 @@ struct dense_output_factory< Stepper , rosenbrock4_dense_output< rosenbrock4_con
     typedef Stepper stepper_type;
     typedef rosenbrock4_controller< stepper_type > controller_type;
     typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::time_type time_type;
     typedef rosenbrock4_dense_output< controller_type > dense_output_type;
 
     dense_output_type operator()( value_type abs_error , value_type rel_error , const stepper_type &stepper )
     {
         return dense_output_type( controller_type( abs_error , rel_error , stepper ) );
     }
+
+    dense_output_type operator()( value_type abs_error , value_type rel_error ,
+                                  time_type max_dt, const stepper_type &stepper )
+    {
+        return dense_output_type( controller_type( abs_error , rel_error , max_dt , stepper ) );
+    }
 };
 
 

include/boost/numeric/odeint/stepper/generation/make_controlled.hpp

@@ -43,6 +43,15 @@ struct controller_factory
     {
         return Controller( abs_error , rel_error , stepper );
     }
+
+    Controller operator()(
+            typename Stepper::value_type abs_error ,
+            typename Stepper::value_type rel_error ,
+            typename Stepper::time_type max_dt ,
+            const Stepper &stepper )
+    {
+        return Controller( abs_error , rel_error , max_dt, stepper );
+    }
 };
 
 
@@ -72,6 +81,19 @@ typename result_of::make_controlled< Stepper >::type make_controlled(
 }
 
 
+template< class Stepper >
+typename result_of::make_controlled< Stepper >::type make_controlled(
+        typename Stepper::value_type abs_error ,
+        typename Stepper::value_type rel_error ,
+        typename Stepper::time_type max_dt ,
+        const Stepper & stepper = Stepper() )
+{
+    typedef Stepper stepper_type;
+    typedef typename result_of::make_controlled< stepper_type >::type controller_type;
+    typedef controller_factory< stepper_type , controller_type > factory_type;
+    factory_type factory;
+    return factory( abs_error , rel_error , max_dt, stepper );
+}
 
 } // odeint
 } // numeric

include/boost/numeric/odeint/stepper/generation/make_dense_output.hpp

@@ -39,6 +39,15 @@ struct dense_output_factory
     {
         return DenseOutput( abs_error , rel_error , stepper );
     }
+
+    DenseOutput operator()(
+            typename Stepper::value_type abs_error ,
+            typename Stepper::value_type rel_error ,
+            typename Stepper::time_type max_dt ,
+            const Stepper &stepper )
+    {
+        return DenseOutput( abs_error , rel_error , max_dt , stepper );
+    }
 };
 
 
@@ -68,6 +77,19 @@ typename result_of::make_dense_output< Stepper >::type make_dense_output(
 }
 
 
+template< class Stepper >
+typename result_of::make_dense_output< Stepper >::type make_dense_output(
+        typename Stepper::value_type abs_error ,
+        typename Stepper::value_type rel_error ,
+        typename Stepper::time_type max_dt ,
+        const Stepper &stepper = Stepper() )
+{
+    typedef Stepper stepper_type;
+    typedef typename result_of::make_dense_output< stepper_type >::type dense_output_type;
+    typedef dense_output_factory< stepper_type , dense_output_type > factory_type;
+    factory_type factory;
+    return factory( abs_error , rel_error , max_dt, stepper );
+}
 
 
 } // odeint

include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp

@@ -27,6 +27,7 @@
 
 #include <boost/numeric/odeint/util/copy.hpp>
 #include <boost/numeric/odeint/util/is_resizeable.hpp>
+#include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
 
 #include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
 
@@ -55,12 +56,20 @@ public:
     typedef rosenbrock4_controller< Stepper > controller_type;
 
 
-    rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 , const stepper_type &stepper = stepper_type() )
-    : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) ,
-      m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
-      m_last_rejected( false )
+    rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 ,
+                            const stepper_type &stepper = stepper_type() )
+        : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) ,
+          m_max_dt( static_cast<time_type>(0) ) ,
+          m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
+          m_last_rejected( false )
     { }
 
+    rosenbrock4_controller( value_type atol, value_type rtol, time_type max_dt,
+                            const stepper_type &stepper = stepper_type() )
+            : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) , m_max_dt( max_dt ) ,
+              m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
+              m_last_rejected( false )
+    { }
 
     value_type error( const state_type &x , const state_type &xold , const state_type &xerr )
     {
@@ -104,6 +113,14 @@ public:
     boost::numeric::odeint::controlled_step_result
     try_step( System sys , const state_type &x , time_type &t , state_type &xout , time_type &dt )
     {
+        if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
+        {
+            // given step size is bigger then max_dt
+            // set limit and return fail
+            dt = m_max_dt;
+            return fail;
+        }
+
         BOOST_USING_STD_MIN();
         BOOST_USING_STD_MAX();
         using std::pow;
@@ -142,7 +159,11 @@ public:
                 min BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) :
                 max BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) );
             t += dt;
-            dt = dt_new;
+            // limit step size to max_dt
+            if( m_max_dt != static_cast<time_type>(0) )
+            {
+                dt = detail::min_abs(m_max_dt, dt_new);
+            }
             m_last_rejected = false;
             return success;
         }
@@ -198,6 +219,7 @@ private:
     wrapped_state_type m_xerr;
     wrapped_state_type m_xnew;
     value_type m_atol , m_rtol;
+    time_type m_max_dt;
     bool m_first_step;
     value_type m_err_old , m_dt_old;
     bool m_last_rejected;

include/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp

@@ -7,7 +7,7 @@
  [end_description]
 
  Copyright 2011-2012 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
 
  Distributed under the Boost Software License, Version 1.0.
@@ -27,6 +27,8 @@
 #include <boost/numeric/odeint/stepper/rosenbrock4_controller.hpp>
 #include <boost/numeric/odeint/util/is_resizeable.hpp>
 
+#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
+
 
 namespace boost {
 namespace numeric {
@@ -73,16 +75,13 @@ public:
     template< class System >
     std::pair< time_type , time_type > do_step( System system )
     {
-        const size_t max_count = 1000;
-
+        failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
         controlled_step_result res = fail;
         m_t_old = m_t;
-        size_t count = 0;
         do
         {
             res = m_stepper.try_step( system , get_current_state() , m_t , get_old_state() , m_dt );
-            if( count++ == max_count )
-                throw std::overflow_error( "rosenbrock4 : too much iterations!");
+            fail_checker();  // check for overflow of failed steps
         }
         while( res == fail );
         m_stepper.stepper().prepare_dense_output();

include/boost/numeric/odeint/util/detail/less_with_sign.hpp

@@ -6,7 +6,7 @@
  Helper function to compare times taking into account the sign of dt
  [end_description]
 
- Copyright 2012-2013 Mario Mulansky
+ Copyright 2012-2015 Mario Mulansky
  Copyright 2012 Karsten Ahnert
 
  Distributed under the Boost Software License, Version 1.0.
@@ -57,7 +57,7 @@ T min_abs( T t1 , T t2 )
 {
     BOOST_USING_STD_MIN();
     BOOST_USING_STD_MAX();
-    if( t1>0 )
+    if( get_unit_value(t1)>0 )
         return min BOOST_PREVENT_MACRO_SUBSTITUTION ( t1 , t2 );
     else
         return max BOOST_PREVENT_MACRO_SUBSTITUTION ( t1 , t2 );
@@ -68,7 +68,7 @@ T max_abs( T t1 , T t2 )
 {
     BOOST_USING_STD_MIN();
     BOOST_USING_STD_MAX();
-    if( t1>0 )
+    if( get_unit_value(t1)>0 )
         return max BOOST_PREVENT_MACRO_SUBSTITUTION ( t1 , t2 );
     else
         return min BOOST_PREVENT_MACRO_SUBSTITUTION ( t1 , t2 );

include/boost/numeric/odeint/util/odeint_error.hpp

@@ -0,0 +1,77 @@
+/*
+ [auto_generated]
+ boost/numeric/odeint/util/odeint_error.hpp
+
+ [begin_description]
+ Runtime Exceptions thrown by odeint
+ [end_description]
+
+ Copyright 2015 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+
+#ifndef BOOST_NUMERIC_ODEINT_UTIL_ODEINT_ERROR_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_UTIL_ODEINT_ERROR_HPP_INCLUDED
+
+#include <stdexcept>
+#include <string>
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+
+/**
+ * \brief Runtime error thrown by odeint
+ */
+class odeint_error : public std::runtime_error
+{
+public:
+    odeint_error(const std::string &s)
+            : std::runtime_error(s)
+    { }
+};
+
+
+/**
+ * \brief Runtime error thrown from integrate routines
+ *
+ * This Error occures when too many iterations are performed in between two
+ * observer calls in the integrate routines.
+ */
+class no_progress_error : public odeint_error
+{
+public:
+    no_progress_error(const std::string &s)
+            : odeint_error(s)
+    { }
+};
+
+
+/**
+ * \brief Runtime error thrown during stepsize adjustment
+ *
+ * This Error occures when too many iterations are performed without finding
+ * an appropriate new step size. This usually indicates non-continuous points
+ * in the ODE.
+ */
+class step_adjustment_error : public odeint_error
+{
+public:
+    step_adjustment_error(const std::string &s)
+            : odeint_error(s)
+    { }
+};
+
+}
+}
+}
+
+
+
+#endif // BOOST_NUMERIC_ODEINT_UTIL_ODEINT_ERROR_HPP_INCLUDED

test/Jamfile.v2

@@ -78,6 +78,8 @@ test-suite "odeint"
     [ run n_step_time_iterator.cpp ]
     [ run times_iterator.cpp ]
     [ run times_time_iterator.cpp ]
+    [ run step_size_limitation.cpp ]
+    [ run integrate_overflow.cpp ]
     [ compile unwrap_boost_reference.cpp ]
     [ compile unwrap_reference.cpp : <cxxflags>-std=c++0x : unwrap_reference_C++11 ]
     [ compile-fail unwrap_reference.cpp : <cxxflags>-std=c++98 : unwrap_reference_C++98 ]

test/integrate.cpp

@@ -97,17 +97,20 @@ struct perform_integrate_const_test
 {
     void operator()( const value_type t_end , const value_type dt )
     {
-        std::cout << "Testing integrate_const with " << typeid( Stepper ).name() << std::endl;
+        // std::cout << "Testing integrate_const with " << typeid( Stepper ).name() << std::endl;
 
         state_type x( 3 , 10.0 ) , x_end( 3 );
 
         std::vector< value_type > times;
 
         size_t steps = integrate_const( Stepper() , lorenz , x , 0.0 , t_end ,
-                         dt , push_back_time( times , x_end ) );
+                                        dt , push_back_time( times , x_end ) );
 
+        std::cout.precision(16);
         std::cout << t_end << " (" << dt << "), " << steps << " , " << times.size() << " , " << 10.0+dt*steps << "=" << x_end[0] << std::endl;
 
+        std::cout << static_cast<int>(floor(t_end/dt)) << " , " << t_end/dt << std::endl;
+
         BOOST_CHECK_EQUAL( static_cast<int>(times.size()) , static_cast<int>(floor(t_end/dt))+1 );
 
         for( size_t i=0 ; i<times.size() ; ++i )
@@ -129,7 +132,7 @@ struct perform_integrate_adaptive_test
 {
     void operator()( const value_type t_end = 10.0 , const value_type dt = 0.03 )
     {
-        std::cout << "Testing integrate_adaptive with " << typeid( Stepper ).name() << std::endl;
+        // std::cout << "Testing integrate_adaptive with " << typeid( Stepper ).name() << std::endl;
 
         state_type x( 3 , 10.0 ) , x_end( 3 );
 
@@ -138,7 +141,7 @@ struct perform_integrate_adaptive_test
         size_t steps = integrate_adaptive( Stepper() , lorenz , x , 0.0 , t_end ,
                                         dt , push_back_time( times , x_end ) );
 
-        std::cout << t_end << " , " << steps << " , " << times.size() << " , " << dt << " , " << 10.0+t_end << "=" << x_end[0] << std::endl;
+        // std::cout << t_end << " , " << steps << " , " << times.size() << " , " << dt << " , " << 10.0+t_end << "=" << x_end[0] << std::endl;
 
         BOOST_CHECK_EQUAL( times.size() , steps+1 );
 
@@ -158,7 +161,7 @@ struct perform_integrate_times_test
 {
     void operator()( const int n = 10 , const int dn=1 , const value_type dt = 0.03 )
     {
-        std::cout << "Testing integrate_times with " << typeid( Stepper ).name() << std::endl;
+        // std::cout << "Testing integrate_times with " << typeid( Stepper ).name() << std::endl;
 
         state_type x( 3 ) , x_end( 3 );
         x[0] = x[1] = x[2] = 10.0;
@@ -194,7 +197,8 @@ struct perform_integrate_n_steps_test
 {
     void operator()( const int n = 200 , const value_type dt = 0.01 )
     {
-        std::cout << "Testing integrate_n_steps with " << typeid( Stepper ).name() << std::endl;
+        // std::cout << "Testing integrate_n_steps with " << typeid( Stepper ).name() << ". ";
+        // std::cout << "dt=" << dt << std::endl;
 
         state_type x( 3 ) , x_end( 3 );
         x[0] = x[1] = x[2] = 10.0;
@@ -208,9 +212,10 @@ struct perform_integrate_n_steps_test
         BOOST_CHECK_EQUAL( static_cast<int>(times.size()) , n+1 );
 
         for( size_t i=0 ; i<times.size() ; ++i )
+        {
             // check if observer was called at times 0,1,2,...
-            BOOST_CHECK_SMALL( times[i] - static_cast< value_type >(i)*dt , 2E-16 );
-
+            BOOST_CHECK_SMALL(times[i] - static_cast< value_type >(i) * dt, 2E-16);
+        }
         // check first, trivial, component
         BOOST_CHECK_SMALL( (10.0 + end_time) - x_end[0] , 1E-6 ); // precision of steppers: 1E-6
 //        BOOST_CHECK_EQUAL( x[1] , x_end[1] );

test/integrate_overflow.cpp

@@ -0,0 +1,200 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/test/integrate_overflow.cpp
+
+ [begin_description]
+ This file tests the overflow exception of the integrate functions.
+ [end_description]
+
+ Copyright 2015 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+
+#define BOOST_TEST_MODULE odeint_integrate_overflow
+
+#include <boost/test/unit_test.hpp>
+
+#include <utility>
+#include <iostream>
+#include <vector>
+
+
+#include <vector>
+#include <cmath>
+#include <iostream>
+
+#include <boost/iterator/counting_iterator.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+
+using namespace boost::numeric::odeint;
+
+
+typedef double value_type;
+typedef std::vector< value_type > state_type;
+
+// the famous lorenz system as usual
+void lorenz( const state_type &x , state_type &dxdt , const value_type t )
+{
+    //const value_type sigma( 10.0 );
+    const value_type R( 28.0 );
+    const value_type b( value_type( 8.0 ) / value_type( 3.0 ) );
+
+    // first component trivial
+    dxdt[0] = 1.0; //sigma * ( x[1] - x[0] );
+    dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+    dxdt[2] = -b * x[2] + x[0] * x[1];
+}
+
+struct push_back_time
+{
+    std::vector< double >& m_times;
+
+    push_back_time( std::vector< double > &times )
+            :  m_times( times ) { }
+
+    void operator()( const state_type &x , double t )
+    {
+        m_times.push_back( t );
+    }
+};
+
+typedef runge_kutta_dopri5<state_type> stepper_type;
+typedef runge_kutta_cash_karp54<state_type> cash_karp_type;
+
+
+BOOST_AUTO_TEST_SUITE( integrate_overflow )
+
+BOOST_AUTO_TEST_CASE( test_integrate_const )
+{
+    state_type x(3, 5.0);
+    std::vector<double> times;
+
+    // check the function signatures with normal stepper
+    integrate_const(stepper_type(), lorenz, x, 0.0, 10.0, 1.0);
+    integrate_const(stepper_type(), lorenz, x, 0.0, 10.0, 1.0, null_observer());
+    // no exceptions expected for normal steppers
+    integrate_const(stepper_type(), lorenz, x, 0.0, 10.0, 1.0, null_observer(), max_step_checker(10));
+
+
+    // check exceptions for controlled stepper
+    integrate_const(make_controlled<stepper_type>(1E-5, 1E-5), lorenz, x, 0.0, 10.0, 1.0);
+    // very small error terms -> standard overflow threshold of 500 should fire an exception
+    BOOST_CHECK_THROW(integrate_const(make_controlled<stepper_type>(1E-15, 1E-15), lorenz, x,
+                                      0.0, 10.0, 1.0, push_back_time(times), max_step_checker()),
+                      std::runtime_error);
+    // small threshold of 10 -> larger error still gives an exception
+    BOOST_CHECK_THROW(integrate_const(make_controlled<stepper_type>(1E-5, 1E-5), lorenz, x,
+                                      0.0, 10.0, 1.0, push_back_time(times), max_step_checker(10)),
+                      std::runtime_error);
+
+    // check exceptions for dense output stepper
+    integrate_const(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x, 0.0, 10.0, 1.0);
+    integrate_const(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x, 0.0, 10.0, 1.0);
+    // very small error terms -> standard overflow threshold of 500 should fire an exception
+    BOOST_CHECK_THROW(integrate_const(make_dense_output<stepper_type>(1E-15, 1E-15), lorenz, x,
+                                      0.0, 10.0, 1.0, push_back_time(times), max_step_checker()),
+                      std::runtime_error);
+    // small threshold of 10 -> larger error still gives an exception
+    BOOST_CHECK_THROW(integrate_const(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x,
+                                      0.0, 10.0, 1.0, push_back_time(times), max_step_checker(10)),
+                      std::runtime_error);
+}
+
+BOOST_AUTO_TEST_CASE( test_integrate_n_steps )
+{
+    state_type x(3, 5.0);
+    std::vector<double> times;
+
+    // check the function signatures with normal stepper
+    integrate_n_steps(stepper_type(), lorenz, x, 0.0, 1.0, 10);
+    integrate_n_steps(stepper_type(), lorenz, x, 0.0, 1.0, 10, null_observer());
+    // no exceptions expected for normal steppers
+    integrate_n_steps(stepper_type(), lorenz, x, 0.0, 1.0, 10, null_observer(), max_step_checker(10));
+
+
+    // check exceptions for controlled stepper
+    integrate_n_steps(make_controlled<stepper_type>(1E-5, 1E-5), lorenz, x, 0.0, 1.0, 10);
+    // very small error terms -> standard overflow threshold of 500 should fire an exception
+    BOOST_CHECK_THROW(integrate_n_steps(make_controlled<stepper_type>(1E-15, 1E-15), lorenz, x,
+                                        0.0, 1.0, 10, push_back_time(times), max_step_checker()),
+                      std::runtime_error);
+    // small threshold of 10 -> larger error still gives an exception
+    BOOST_CHECK_THROW(integrate_n_steps(make_controlled<stepper_type>(1E-5, 1E-5), lorenz, x,
+                                        0.0, 1.0, 10, push_back_time(times), max_step_checker(10)),
+                      std::runtime_error);
+
+    // check exceptions for dense output stepper
+    integrate_n_steps(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x, 0.0, 1.0, 10);
+    integrate_n_steps(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x, 0.0, 1.0, 10, push_back_time(times));
+    // very small error terms -> standard overflow threshold of 500 should fire an exception
+    BOOST_CHECK_THROW(integrate_n_steps(make_dense_output<stepper_type>(1E-15, 1E-15), lorenz, x,
+                                        0.0, 1.0, 10, push_back_time(times), max_step_checker()),
+                      std::runtime_error);
+    // small threshold of 10 -> larger error still gives an exception
+    BOOST_CHECK_THROW(integrate_n_steps(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x,
+                                        0.0, 1.0, 10, push_back_time(times), max_step_checker(10)),
+                      std::runtime_error);
+}
+
+BOOST_AUTO_TEST_CASE( test_integrate_times )
+{
+    state_type x(3, 5.0);
+    std::vector<double> times;
+    boost::counting_iterator<int> t0(0);
+    boost::counting_iterator<int> t1(10);
+
+    // check the function signatures with normal stepper
+    integrate_times(stepper_type(), lorenz, x, t0, t1, 1.0 , push_back_time(times));
+    // no exceptions expected for big enough step size
+    integrate_times(stepper_type(), lorenz, x, t0, t1, 1.0 , push_back_time(times), max_step_checker(10));
+    // if dt*max_steps < observer time difference we expect an exception
+    BOOST_CHECK_THROW(integrate_times(stepper_type(), lorenz, x, t0, t1, 0.01, push_back_time(times),
+                                      max_step_checker(10)),
+                      std::runtime_error);
+
+    // check exceptions for controlled stepper
+    // no exception if no checker is provided
+    integrate_times(make_controlled<stepper_type>(1E-5, 1E-5), lorenz, x, t0, t1, 1.0 , push_back_time(times));
+    // very small error terms -> standard overflow threshold of 500 should fire an exception
+    BOOST_CHECK_THROW(integrate_times(make_controlled<stepper_type>(1E-15, 1E-15), lorenz, x,
+                                      t0, t1, 1.0 , push_back_time(times), max_step_checker()),
+                      std::runtime_error);
+
+    // small threshold of 10 -> larger error still gives an exception
+    BOOST_CHECK_THROW(integrate_times(make_controlled<stepper_type>(1E-5, 1E-5), lorenz, x,
+                                      t0, t1, 1.0 , push_back_time(times), max_step_checker(10)),
+                      std::runtime_error);
+
+    // check cash karp controlled stepper
+    // no exception if no checker is provided
+    integrate_times(make_controlled<cash_karp_type>(1E-5, 1E-5), lorenz, x, t0, t1, 1.0 , push_back_time(times));
+    // very small error terms -> standard overflow threshold of 500 should fire an exception
+    BOOST_CHECK_THROW(integrate_times(make_controlled<cash_karp_type>(1E-15, 1E-15), lorenz, x,
+                      t0, t1, 1.0 , push_back_time(times), max_step_checker()),
+    std::runtime_error);
+
+    // small threshold of 10 -> larger error still gives an exception
+    BOOST_CHECK_THROW(integrate_times(make_controlled<cash_karp_type>(1E-5, 1E-5), lorenz, x,
+                      t0, t1, 1.0 , push_back_time(times), max_step_checker(10)),
+    std::runtime_error);
+
+    // check exceptions for dense output stepper
+    integrate_times(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x, t0, t1, 1.0 ,
+                    push_back_time(times));
+    // very small error terms -> standard overflow threshold of 500 should fire an exception
+    BOOST_CHECK_THROW(integrate_times(make_dense_output<stepper_type>(1E-15, 1E-15), lorenz, x,
+                                      t0, t1, 1.0 , push_back_time(times), max_step_checker()),
+                      std::runtime_error);
+    // small threshold of 10 -> larger error still gives an exception
+    BOOST_CHECK_THROW(integrate_times(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x,
+                                      t0, t1, 1.0 , push_back_time(times), max_step_checker(10)),
+                      std::runtime_error);
+}
+
+BOOST_AUTO_TEST_SUITE_END()

test/step_size_limitation.cpp

@@ -0,0 +1,268 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/test/step_size_limitation.cpp
+
+ [begin_description]
+ Tests the step size limitation functionality
+ [end_description]
+
+ Copyright 2015 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#define BOOST_TEST_MODULE odeint_integrate_times
+
+#include <boost/test/unit_test.hpp>
+
+#include <utility>
+#include <iostream>
+#include <vector>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+typedef double value_type;
+typedef std::vector< value_type > state_type;
+
+/***********************************************
+ * first part of the tests: explicit methods
+ ***********************************************
+ */
+
+void damped_osc( const state_type &x , state_type &dxdt , const value_type t )
+{
+    const value_type gam( 0.1);
+
+    dxdt[0] = x[1];
+    dxdt[1] = -x[0] - gam*x[1];
+}
+
+
+struct push_back_time
+{
+    std::vector< double >& m_times;
+
+    push_back_time( std::vector< double > &times )
+    :  m_times( times ) { }
+
+    template<typename State>
+    void operator()( const State &x , double t )
+    {
+        m_times.push_back( t );
+    }
+};
+
+BOOST_AUTO_TEST_SUITE( step_size_limitation_test )
+
+BOOST_AUTO_TEST_CASE( test_step_adjuster )
+{
+    // first use adjuster without step size limitation
+    default_step_adjuster<double, double> step_adjuster;
+    const double dt = 0.1;
+    double dt_new = step_adjuster.decrease_step(dt, 1.5, 2);
+    BOOST_CHECK(dt_new < dt*2.0/3.0);
+
+    dt_new = step_adjuster.increase_step(dt, 0.8, 1);
+    // for errors > 0.5 no increase is performed
+    BOOST_CHECK(dt_new == dt);
+
+    dt_new = step_adjuster.increase_step(dt, 0.4, 1);
+    // smaller errors should lead to step size increase
+    std::cout << dt_new << std::endl;
+    BOOST_CHECK(dt_new > dt);
+
+
+    // now test with step size limitation max_dt = 0.1
+    default_step_adjuster<double, double>
+        limited_adjuster(dt);
+
+    dt_new = limited_adjuster.decrease_step(dt, 1.5, 2);
+    // decreasing works as before
+    BOOST_CHECK(dt_new < dt*2.0/3.0);
+
+    dt_new = limited_adjuster.decrease_step(2*dt, 1.1, 2);
+    // decreasing a large step size should give max_dt
+    BOOST_CHECK(dt_new == dt);
+
+    dt_new = limited_adjuster.increase_step(dt, 0.8, 1);
+    // for errors > 0.5 no increase is performed, still valid
+    BOOST_CHECK(dt_new == dt);
+
+    dt_new = limited_adjuster.increase_step(dt, 0.4, 1);
+    // but even for smaller errors, we should at most get 0.1
+    BOOST_CHECK_EQUAL(dt_new, dt);
+
+    dt_new = limited_adjuster.increase_step(0.9*dt, 0.1, 1);
+    std::cout << dt_new << std::endl;
+    // check that we don't increase beyond max_dt
+    BOOST_CHECK(dt_new == dt);
+}
+
+
+template<class Stepper>
+void test_explicit_stepper(Stepper stepper, const double max_dt)
+{
+    state_type x( 2 );
+    x[0] = x[1] = 10.0;
+    const size_t steps = 100;
+
+    std::vector<double> times;
+
+    integrate_adaptive(stepper, damped_osc, x, 0.0, steps*max_dt, max_dt, push_back_time(times));
+
+    BOOST_CHECK_EQUAL(times.size(), steps+1);
+    // check that dt remains at exactly max_dt
+    for( size_t i=0 ; i<times.size() ; ++i )
+        // check if observer was called at times 0,1,2,...
+        BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15);
+    times.clear();
+
+    // this should also work when we provide some bigger initial dt
+    integrate_adaptive(stepper, damped_osc, x, 0.0, steps*max_dt, 10*max_dt, push_back_time(times));
+
+    BOOST_CHECK_EQUAL(times.size(), steps+1);
+    // check that dt remains at exactly max_dt
+    for( size_t i=0 ; i<times.size() ; ++i )
+        // check if observer was called at times 0,1,2,...
+        BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15);
+    times.clear();
+}
+
+
+BOOST_AUTO_TEST_CASE( test_controlled )
+{
+    const double max_dt = 0.01;
+
+    test_explicit_stepper(make_controlled(1E-2, 1E-2, max_dt,
+                                            runge_kutta_dopri5<state_type>()),
+                            max_dt);
+    test_explicit_stepper(make_controlled(1E-2, 1E-2, -max_dt,
+                                            runge_kutta_dopri5<state_type>()),
+                            -max_dt);
+
+    test_explicit_stepper(make_controlled(1E-2, 1E-2, max_dt,
+                                            runge_kutta_cash_karp54<state_type>()),
+                            max_dt);
+    test_explicit_stepper(make_controlled(1E-2, 1E-2, -max_dt,
+                                            runge_kutta_cash_karp54<state_type>()),
+                            -max_dt);
+
+    test_explicit_stepper(bulirsch_stoer<state_type>(1E-2, 1E-2, 1.0, 1.0, max_dt),
+                            max_dt);
+    test_explicit_stepper(bulirsch_stoer<state_type>(1E-2, 1E-2, 1.0, 1.0, -max_dt),
+                            -max_dt);
+}
+
+
+BOOST_AUTO_TEST_CASE( test_dense_out )
+{
+    const double max_dt = 0.01;
+    test_explicit_stepper(make_dense_output(1E-2, 1E-2, max_dt,
+                                             runge_kutta_dopri5<state_type>()),
+                          max_dt);
+    test_explicit_stepper(make_dense_output(1E-2, 1E-2, -max_dt,
+                                            runge_kutta_dopri5<state_type>()),
+                          -max_dt);
+
+    test_explicit_stepper(bulirsch_stoer_dense_out<state_type>(1E-2, 1E-2, 1, 1, max_dt),
+                          max_dt);
+
+    test_explicit_stepper(bulirsch_stoer_dense_out<state_type>(1E-2, 1E-2, 1, 1, -max_dt),
+                          -max_dt);
+}
+
+
+/***********************************************
+ * second part of the tests: implicit Rosenbrock
+ ***********************************************
+ */
+
+typedef boost::numeric::ublas::vector< value_type > vector_type;
+typedef boost::numeric::ublas::matrix< value_type > matrix_type;
+
+
+// harmonic oscillator, analytic solution x[0] = sin( t )
+struct osc_rhs
+{
+    void operator()( const vector_type &x , vector_type &dxdt , const value_type &t ) const
+    {
+        dxdt( 0 ) = x( 1 );
+        dxdt( 1 ) = -x( 0 );
+    }
+};
+
+struct osc_jacobi
+{
+    void operator()( const vector_type &x , matrix_type &jacobi , const value_type &t , vector_type &dfdt ) const
+    {
+        jacobi( 0 , 0 ) = 0;
+        jacobi( 0 , 1 ) = 1;
+        jacobi( 1 , 0 ) = -1;
+        jacobi( 1 , 1 ) = 0;
+        dfdt( 0 ) = 0.0;
+        dfdt( 1 ) = 0.0;
+    }
+};
+
+
+template<class Stepper>
+void test_rosenbrock_stepper(Stepper stepper, const double max_dt)
+{
+    vector_type x( 2 );
+    x(0) = x(1) = 10.0;
+    const size_t steps = 100;
+
+    std::vector<double> times;
+
+    integrate_adaptive(stepper,
+                       std::make_pair(osc_rhs(), osc_jacobi()),
+                       x, 0.0, steps*max_dt, max_dt, push_back_time(times));
+
+    BOOST_CHECK_EQUAL(times.size(), steps+1);
+    // check that dt remains at exactly max_dt
+    for( size_t i=0 ; i<times.size() ; ++i )
+        // check if observer was called at times 0,1,2,...
+        BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15);
+    times.clear();
+
+    // this should also work when we provide some bigger initial dt
+    integrate_adaptive(stepper,
+                       std::make_pair(osc_rhs(), osc_jacobi()),
+                       x, 0.0, steps*max_dt, 10*max_dt, push_back_time(times));
+
+    BOOST_CHECK_EQUAL(times.size(), steps+1);
+    // check that dt remains at exactly max_dt
+    for( size_t i=0 ; i<times.size() ; ++i )
+        // check if observer was called at times 0,1,2,...
+        BOOST_CHECK_SMALL( times[i] - static_cast<double>(i)*max_dt , 1E-15);
+    times.clear();
+}
+
+
+BOOST_AUTO_TEST_CASE( test_controlled_rosenbrock )
+{
+    const double max_dt = 0.01;
+
+    test_rosenbrock_stepper(make_controlled(1E-2, 1E-2, max_dt, rosenbrock4<value_type>()),
+                            max_dt);
+    test_rosenbrock_stepper(make_controlled(1E-2, 1E-2, -max_dt, rosenbrock4<value_type>()),
+                            -max_dt);
+}
+
+
+BOOST_AUTO_TEST_CASE( test_dense_out_rosenbrock )
+{
+    const double max_dt = 0.01;
+
+    test_rosenbrock_stepper(make_dense_output(1E-2, 1E-2, max_dt, rosenbrock4<value_type>()),
+                            max_dt);
+    test_rosenbrock_stepper(make_dense_output(1E-2, 1E-2, -max_dt, rosenbrock4<value_type>()),
+                            -max_dt);
+}
+
+BOOST_AUTO_TEST_SUITE_END()