diff options
Diffstat (limited to 'boost/numeric/odeint/stepper')
11 files changed, 319 insertions, 118 deletions
diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp index e71008a993..4b908333ba 100644 --- a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp +++ b/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; }; diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp index 3018524ec9..d876ca3d36 100644 --- a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp +++ b/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; }; diff --git a/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp b/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp index 509192482c..8ae627fe1c 100644 --- a/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp +++ b/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,59 +411,35 @@ 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 diff --git a/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp b/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp index b685d92bcf..94abc5af99 100644 --- a/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp +++ b/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(); diff --git a/boost/numeric/odeint/stepper/generation/generation_controlled_runge_kutta.hpp b/boost/numeric/odeint/stepper/generation/generation_controlled_runge_kutta.hpp index 54b7746c6e..bad2bedfef 100644 --- a/boost/numeric/odeint/stepper/generation/generation_controlled_runge_kutta.hpp +++ b/boost/numeric/odeint/stepper/generation/generation_controlled_runge_kutta.hpp @@ -34,15 +34,23 @@ 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 diff --git a/boost/numeric/odeint/stepper/generation/generation_dense_output_runge_kutta.hpp b/boost/numeric/odeint/stepper/generation/generation_dense_output_runge_kutta.hpp index 214d09c1a4..276358a12e 100644 --- a/boost/numeric/odeint/stepper/generation/generation_dense_output_runge_kutta.hpp +++ b/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 ) ); } }; diff --git a/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp b/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp index d5322c8d73..366cb5e9c4 100644 --- a/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp +++ b/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 ) ); + } }; diff --git a/boost/numeric/odeint/stepper/generation/make_controlled.hpp b/boost/numeric/odeint/stepper/generation/make_controlled.hpp index 603978fec8..61bc8f1920 100644 --- a/boost/numeric/odeint/stepper/generation/make_controlled.hpp +++ b/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 diff --git a/boost/numeric/odeint/stepper/generation/make_dense_output.hpp b/boost/numeric/odeint/stepper/generation/make_dense_output.hpp index a28e31b5d2..fff3590525 100644 --- a/boost/numeric/odeint/stepper/generation/make_dense_output.hpp +++ b/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 diff --git a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp index df4e6c48b6..0e5edd32c0 100644 --- a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp +++ b/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; diff --git a/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp b/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp index 8a866c562e..6695ba6a97 100644 --- a/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp +++ b/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(); |