add max_dt support to generation functions + test

generate function now support additional max_dt parameter for setting the set
size limit.
Added a test case to check limiter behavior for controlled and dense out
integration.

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

@@ -93,17 +93,21 @@ public:
     }
 
 
-    time_type decrease_step(const time_type dt, const value_type error, const int error_order) const
+    time_type decrease_step(time_type dt, const value_type error, const int error_order) const
     {
         BOOST_USING_STD_MIN();
         BOOST_USING_STD_MAX();
         using std::pow;
 
-        return dt * max
+        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
@@ -113,8 +117,7 @@ public:
         using std::pow;
 
         // adjust the size if dt is smaller than max_dt (providede max_dt is not zero)
-        if(error < 0.5 && (m_max_dt == static_cast<time_type >(0) ||
-                           detail::less_with_sign(dt, m_max_dt, dt)))
+        if(error < 0.5)
         {
             // error should be > 0
             error = max BOOST_PREVENT_MACRO_SUBSTITUTION (
@@ -124,12 +127,22 @@ public:
             //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);
-            // limit to maximal stepsize
-            dt = detail::max_abs(dt, m_max_dt);
+            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:
 
@@ -380,6 +393,12 @@ 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_error_checker.check_step_size_limit(dt) )
+        {
+            // given dt was above step size limit - adjust and return fail;
+            dt = m_error_checker.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 ) );
 
@@ -710,10 +729,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_error_checker.check_step_size_limit(dt) )
+        {
+            // given dt was above step size limit - adjust and return fail;
+            dt = m_error_checker.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 ) );
 

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

@@ -35,16 +35,21 @@ struct controller_factory< Stepper , controlled_runge_kutta< Stepper > >
     typedef controlled_runge_kutta< stepper_type > controller_type;
     typedef typename controller_type::error_checker_type error_checker_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 );
     }
+
+    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 , 1, 1, max_dt) , stepper );
+    }
 };
 
 
-
-
 } // odeint
 } // numeric
 } // boost

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

@@ -34,12 +34,21 @@ struct dense_output_factory< Stepper , dense_output_runge_kutta< controlled_rung
     typedef controlled_runge_kutta< stepper_type > controller_type;
     typedef typename controller_type::error_checker_type error_checker_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 ) );
     }
+
+    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 , 1, 1, 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/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.

test/Jamfile.v2

@@ -78,6 +78,7 @@ test-suite "odeint"
     [ run n_step_time_iterator.cpp ]
     [ run times_iterator.cpp ]
     [ run times_time_iterator.cpp ]
+    [ run step_size_limitation.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/step_size_limitation.cpp

@@ -0,0 +1,179 @@
+/*
+ [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;
+
+
+void lorenz( const state_type &x , state_type &dxdt , const value_type t )
+{
+    BOOST_CHECK( t >= 0.0 );
+
+    const value_type sigma( 10.0 );
+    const value_type R( 28.0 );
+    const value_type b( value_type( 8.0 ) / value_type( 3.0 ) );
+
+    dxdt[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 );
+    }
+};
+
+BOOST_AUTO_TEST_SUITE( step_size_limitation_test )
+
+BOOST_AUTO_TEST_CASE( test_error_checker )
+{
+    // first use checker without step size limitation
+    default_error_checker<double, double, range_algebra, default_operations> checker;
+    const double dt = 0.1;
+    double dt_new = checker.decrease_step(dt, 1.5, 2);
+    BOOST_CHECK(dt_new < dt*2.0/3.0);
+
+    dt_new = checker.increase_step(dt, 0.8, 1);
+    // for errors > 0.5 no increase is performed
+    BOOST_CHECK(dt_new == dt);
+
+    dt_new = checker.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_error_checker<double, double, range_algebra, default_operations>
+        limited_checker(1E-6, 1E-6, 1, 1, dt);
+
+    dt_new = limited_checker.decrease_step(dt, 1.5, 2);
+    // decreasing works as before
+    BOOST_CHECK(dt_new < dt*2.0/3.0);
+
+    dt_new = limited_checker.decrease_step(2*dt, 1.1, 2);
+    // decreasing a large step size should give max_dt
+    BOOST_CHECK(dt_new == dt);
+
+    dt_new = limited_checker.increase_step(dt, 0.8, 1);
+    // for errors > 0.5 no increase is performed, still valid
+    BOOST_CHECK(dt_new == dt);
+
+    dt_new = limited_checker.increase_step(dt, 0.4, 1);
+    // but even for smaller errors, we should at most get 0.1
+    BOOST_CHECK(dt_new == dt);
+
+    dt_new = limited_checker.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);
+}
+
+
+BOOST_AUTO_TEST_CASE( test_controlled )
+{
+    state_type x( 3 );
+    x[0] = x[1] = x[2] = 10.0;
+    const double max_dt = 0.01;
+
+    std::vector<double> times;
+
+    integrate_adaptive(make_controlled(1E-2, 1E-2, max_dt, runge_kutta_dopri5<state_type>()),
+                       lorenz, x, 0.0, 1.0, max_dt, push_back_time(times));
+    // 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(make_controlled(1E-2, 1E-2, max_dt, runge_kutta_dopri5<state_type>()),
+                       lorenz, x, 0.0, 1.0, 0.1, push_back_time(times));
+    // 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();
+
+    // check this behavior with cash karp
+    integrate_adaptive(make_controlled(1E-2, 1E-2, max_dt, runge_kutta_cash_karp54<state_type>()),
+                       lorenz, x, 0.0, 1.0, max_dt, push_back_time(times));
+    // 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(make_controlled(1E-2, 1E-2, max_dt, runge_kutta_cash_karp54<state_type>()),
+                       lorenz, x, 0.0, 1.0, 0.1, push_back_time(times));
+    // 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_dense_out )
+{
+    state_type x( 3 );
+    x[0] = x[1] = x[2] = 10.0;
+    const double max_dt = 0.01;
+
+    std::vector<double> times;
+
+    integrate_adaptive(make_dense_output(1E-2, 1E-2, max_dt, runge_kutta_dopri5<state_type>()),
+                      lorenz, x, 0.0, 1.0, max_dt, push_back_time(times));
+    // 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(make_dense_output(1E-2, 1E-2, max_dt, runge_kutta_dopri5<state_type>()),
+                       lorenz, x, 0.0, 1.0, 0.1, push_back_time(times));
+    // 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_SUITE_END()