diff options
Diffstat (limited to 'boost/numeric/odeint/integrate')
10 files changed, 669 insertions, 144 deletions
diff --git a/boost/numeric/odeint/integrate/check_adapter.hpp b/boost/numeric/odeint/integrate/check_adapter.hpp new file mode 100644 index 0000000000..3d3ebd6c88 --- /dev/null +++ b/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
\ No newline at end of file diff --git a/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp b/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp index 743e57709c..7516d44087 100644 --- a/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp +++ b/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; } diff --git a/boost/numeric/odeint/integrate/detail/integrate_const.hpp b/boost/numeric/odeint/integrate/detail/integrate_const.hpp index 312564ff16..7a86b32fa6 100644 --- a/boost/numeric/odeint/integrate/detail/integrate_const.hpp +++ b/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 ); diff --git a/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp b/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp index 3c1d171620..2ef490d592 100644 --- a/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp +++ b/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; diff --git a/boost/numeric/odeint/integrate/detail/integrate_times.hpp b/boost/numeric/odeint/integrate/detail/integrate_times.hpp index d5446ba590..2e27990412 100644 --- a/boost/numeric/odeint/integrate/detail/integrate_times.hpp +++ b/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 ) diff --git a/boost/numeric/odeint/integrate/integrate_adaptive.hpp b/boost/numeric/odeint/integrate/integrate_adaptive.hpp index 12bdc8237d..09997142fc 100644 --- a/boost/numeric/odeint/integrate/integrate_adaptive.hpp +++ b/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 diff --git a/boost/numeric/odeint/integrate/integrate_const.hpp b/boost/numeric/odeint/integrate/integrate_const.hpp index ba2d891e97..fae683f978 100644 --- a/boost/numeric/odeint/integrate/integrate_const.hpp +++ b/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 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, 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 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, 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 > +* \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() ); + 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 > +* \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() ); + 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 diff --git a/boost/numeric/odeint/integrate/integrate_n_steps.hpp b/boost/numeric/odeint/integrate/integrate_n_steps.hpp index f85c453bd2..7f3a49bddc 100644 --- a/boost/numeric/odeint/integrate/integrate_n_steps.hpp +++ b/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,35 +34,80 @@ 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_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( + 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; + 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() ); + 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. */ diff --git a/boost/numeric/odeint/integrate/integrate_times.hpp b/boost/numeric/odeint/integrate/integrate_times.hpp index c0ecda6c6d..79fba4f1b3 100644 --- a/boost/numeric/odeint/integrate/integrate_times.hpp +++ b/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 × , 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 × , 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. */ diff --git a/boost/numeric/odeint/integrate/max_step_checker.hpp b/boost/numeric/odeint/integrate/max_step_checker.hpp new file mode 100644 index 0000000000..6808a57bdb --- /dev/null +++ b/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
\ No newline at end of file |