Chris@16: /* [auto_generated] Chris@16: boost/numeric/odeint/stepper/controlled_runge_kutta.hpp Chris@16: Chris@16: [begin_description] Chris@16: The default controlled stepper which can be used with all explicit Runge-Kutta error steppers. Chris@16: [end_description] Chris@16: Chris@101: Copyright 2010-2013 Karsten Ahnert Chris@101: Copyright 2010-2013 Mario Mulansky Chris@101: Copyright 2012 Christoph Koke Chris@16: Chris@16: Distributed under the Boost Software License, Version 1.0. Chris@16: (See accompanying file LICENSE_1_0.txt or Chris@16: copy at http://www.boost.org/LICENSE_1_0.txt) Chris@16: */ Chris@16: Chris@16: Chris@16: #ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED Chris@16: #define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED Chris@16: Chris@16: Chris@16: Chris@16: #include Chris@16: Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: Chris@16: #include Chris@16: #include Chris@101: #include Chris@16: Chris@16: #include Chris@16: #include Chris@16: Chris@16: namespace boost { Chris@16: namespace numeric { Chris@16: namespace odeint { Chris@16: Chris@16: Chris@16: template Chris@16: < Chris@16: class Value , Chris@101: class Algebra , Chris@101: class Operations Chris@16: > Chris@16: class default_error_checker Chris@16: { Chris@16: public: Chris@16: Chris@16: typedef Value value_type; Chris@16: typedef Algebra algebra_type; Chris@16: typedef Operations operations_type; Chris@16: Chris@16: default_error_checker( Chris@16: value_type eps_abs = static_cast< value_type >( 1.0e-6 ) , Chris@16: value_type eps_rel = static_cast< value_type >( 1.0e-6 ) , Chris@16: value_type a_x = static_cast< value_type >( 1 ) , Chris@16: value_type a_dxdt = static_cast< value_type >( 1 ) ) Chris@16: : m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt ) Chris@16: { } Chris@16: Chris@16: Chris@16: template< class State , class Deriv , class Err , class Time > Chris@16: value_type error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const Chris@16: { Chris@16: return error( algebra_type() , x_old , dxdt_old , x_err , dt ); Chris@16: } Chris@16: Chris@16: template< class State , class Deriv , class Err , class Time > Chris@16: value_type error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const Chris@16: { Chris@16: // this overwrites x_err ! Chris@16: algebra.for_each3( x_err , x_old , dxdt_old , Chris@16: typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * get_unit_value( dt ) ) ); Chris@16: Chris@101: // value_type res = algebra.reduce( x_err , Chris@101: // typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0 ) ); Chris@101: return algebra.norm_inf( x_err ); Chris@16: } Chris@16: Chris@16: private: Chris@16: Chris@16: value_type m_eps_abs; Chris@16: value_type m_eps_rel; Chris@16: value_type m_a_x; Chris@16: value_type m_a_dxdt; Chris@16: Chris@16: }; Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: /* Chris@16: * error stepper category dispatcher Chris@16: */ Chris@16: template< Chris@16: class ErrorStepper , Chris@16: class ErrorChecker = default_error_checker< typename ErrorStepper::value_type , Chris@16: typename ErrorStepper::algebra_type , Chris@16: typename ErrorStepper::operations_type > , Chris@16: class Resizer = typename ErrorStepper::resizer_type , Chris@16: class ErrorStepperCategory = typename ErrorStepper::stepper_category Chris@16: > Chris@16: class controlled_runge_kutta ; Chris@16: Chris@16: Chris@16: Chris@16: /* Chris@16: * explicit stepper version Chris@16: * Chris@16: * this class introduces the following try_step overloads Chris@16: * try_step( sys , x , t , dt ) Chris@16: * try_step( sys , x , dxdt , t , dt ) Chris@16: * try_step( sys , in , t , out , dt ) Chris@16: * try_step( sys , in , dxdt , t , out , dt ) Chris@16: */ Chris@16: /** Chris@16: * \brief Implements step size control for Runge-Kutta steppers with error Chris@16: * estimation. Chris@16: * Chris@16: * This class implements the step size control for standard Runge-Kutta Chris@16: * steppers with error estimation. Chris@16: * Chris@16: * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept. Chris@16: * \tparam ErrorChecker The error checker Chris@16: * \tparam Resizer The resizer policy type. Chris@16: */ Chris@16: template< Chris@16: class ErrorStepper , Chris@16: class ErrorChecker , Chris@16: class Resizer Chris@16: > Chris@16: class controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > Chris@16: { Chris@16: Chris@16: public: Chris@16: Chris@16: typedef ErrorStepper stepper_type; Chris@16: typedef typename stepper_type::state_type state_type; Chris@16: typedef typename stepper_type::value_type value_type; Chris@16: typedef typename stepper_type::deriv_type deriv_type; Chris@16: typedef typename stepper_type::time_type time_type; Chris@16: typedef typename stepper_type::algebra_type algebra_type; Chris@16: typedef typename stepper_type::operations_type operations_type; Chris@16: typedef Resizer resizer_type; Chris@16: typedef ErrorChecker error_checker_type; Chris@16: typedef explicit_controlled_stepper_tag stepper_category; Chris@16: Chris@16: #ifndef DOXYGEN_SKIP Chris@16: typedef typename stepper_type::wrapped_state_type wrapped_state_type; Chris@16: typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type; Chris@16: Chris@16: typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > controlled_stepper_type; Chris@16: #endif //DOXYGEN_SKIP Chris@16: Chris@16: Chris@16: /** Chris@16: * \brief Constructs the controlled Runge-Kutta stepper. Chris@16: * \param error_checker An instance of the error checker. Chris@16: * \param stepper An instance of the underlying stepper. Chris@16: */ Chris@16: controlled_runge_kutta( Chris@16: const error_checker_type &error_checker = error_checker_type( ) , Chris@16: const stepper_type &stepper = stepper_type( ) Chris@16: ) Chris@16: : m_stepper( stepper ) , m_error_checker( error_checker ) Chris@16: { } Chris@16: Chris@16: Chris@16: Chris@16: /* Chris@16: * Version 1 : try_step( sys , x , t , dt ) Chris@16: * Chris@16: * The overloads are needed to solve the forwarding problem Chris@16: */ Chris@16: /** Chris@16: * \brief Tries to perform one step. Chris@16: * Chris@16: * This method tries to do one step with step size dt. If the error estimate Chris@16: * is to large, the step is rejected and the method returns fail and the Chris@16: * step size dt is reduced. If the error estimate is acceptably small, the Chris@16: * step is performed, success is returned and dt might be increased to make Chris@16: * the steps as large as possible. This method also updates t if a step is Chris@16: * performed. Chris@16: * Chris@16: * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Chris@16: * Simple System concept. Chris@16: * \param x The state of the ODE which should be solved. Overwritten if Chris@16: * the step is successful. Chris@16: * \param t The value of the time. Updated if the step is successful. Chris@16: * \param dt The step size. Updated. Chris@16: * \return success if the step was accepted, fail otherwise. Chris@16: */ Chris@16: template< class System , class StateInOut > Chris@16: controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt ) Chris@16: { Chris@16: return try_step_v1( system , x , t, dt ); Chris@16: } Chris@16: Chris@16: /** Chris@16: * \brief Tries to perform one step. Solves the forwarding problem and Chris@16: * allows for using boost range as state_type. Chris@16: * Chris@16: * This method tries to do one step with step size dt. If the error estimate Chris@16: * is to large, the step is rejected and the method returns fail and the Chris@16: * step size dt is reduced. If the error estimate is acceptably small, the Chris@16: * step is performed, success is returned and dt might be increased to make Chris@16: * the steps as large as possible. This method also updates t if a step is Chris@16: * performed. Chris@16: * Chris@16: * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Chris@16: * Simple System concept. Chris@16: * \param x The state of the ODE which should be solved. Overwritten if Chris@16: * the step is successful. Can be a boost range. Chris@16: * \param t The value of the time. Updated if the step is successful. Chris@16: * \param dt The step size. Updated. Chris@16: * \return success if the step was accepted, fail otherwise. Chris@16: */ Chris@16: template< class System , class StateInOut > Chris@16: controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt ) Chris@16: { Chris@16: return try_step_v1( system , x , t, dt ); Chris@16: } Chris@16: Chris@16: Chris@16: Chris@16: /* Chris@16: * Version 2 : try_step( sys , x , dxdt , t , dt ) Chris@16: * Chris@16: * this version does not solve the forwarding problem, boost.range can not be used Chris@16: */ Chris@16: /** Chris@16: * \brief Tries to perform one step. Chris@16: * Chris@16: * This method tries to do one step with step size dt. If the error estimate Chris@16: * is to large, the step is rejected and the method returns fail and the Chris@16: * step size dt is reduced. If the error estimate is acceptably small, the Chris@16: * step is performed, success is returned and dt might be increased to make Chris@16: * the steps as large as possible. This method also updates t if a step is Chris@16: * performed. Chris@16: * Chris@16: * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Chris@16: * Simple System concept. Chris@16: * \param x The state of the ODE which should be solved. Overwritten if Chris@16: * the step is successful. Chris@16: * \param dxdt The derivative of state. Chris@16: * \param t The value of the time. Updated if the step is successful. Chris@16: * \param dt The step size. Updated. Chris@16: * \return success if the step was accepted, fail otherwise. Chris@16: */ Chris@16: template< class System , class StateInOut , class DerivIn > Chris@16: controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt ) Chris@16: { Chris@16: m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); Chris@16: controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt ); Chris@16: if( res == success ) Chris@16: { Chris@16: boost::numeric::odeint::copy( m_xnew.m_v , x ); Chris@16: } Chris@16: return res; Chris@16: } Chris@16: Chris@16: /* Chris@16: * Version 3 : try_step( sys , in , t , out , dt ) Chris@16: * Chris@16: * this version does not solve the forwarding problem, boost.range can not be used Chris@16: * Chris@16: * the disable is needed to avoid ambiguous overloads if state_type = time_type Chris@16: */ Chris@16: /** Chris@16: * \brief Tries to perform one step. Chris@16: * Chris@16: * \note This method is disabled if state_type=time_type to avoid ambiguity. Chris@16: * Chris@16: * This method tries to do one step with step size dt. If the error estimate Chris@16: * is to large, the step is rejected and the method returns fail and the Chris@16: * step size dt is reduced. If the error estimate is acceptably small, the Chris@16: * step is performed, success is returned and dt might be increased to make Chris@16: * the steps as large as possible. This method also updates t if a step is Chris@16: * performed. Chris@16: * Chris@16: * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Chris@16: * Simple System concept. Chris@16: * \param in The state of the ODE which should be solved. Chris@16: * \param t The value of the time. Updated if the step is successful. Chris@16: * \param out Used to store the result of the step. Chris@16: * \param dt The step size. Updated. Chris@16: * \return success if the step was accepted, fail otherwise. Chris@16: */ Chris@16: template< class System , class StateIn , class StateOut > Chris@16: typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type Chris@16: try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt ) Chris@16: { Chris@16: typename odeint::unwrap_reference< System >::type &sys = system; Chris@16: m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ); Chris@16: sys( in , m_dxdt.m_v , t ); Chris@16: return try_step( system , in , m_dxdt.m_v , t , out , dt ); Chris@16: } Chris@16: Chris@16: Chris@16: /* Chris@16: * Version 4 : try_step( sys , in , dxdt , t , out , dt ) Chris@16: * Chris@16: * this version does not solve the forwarding problem, boost.range can not be used Chris@16: */ Chris@16: /** Chris@16: * \brief Tries to perform one step. Chris@16: * Chris@16: * This method tries to do one step with step size dt. If the error estimate Chris@16: * is to large, the step is rejected and the method returns fail and the Chris@16: * step size dt is reduced. If the error estimate is acceptably small, the Chris@16: * step is performed, success is returned and dt might be increased to make Chris@16: * the steps as large as possible. This method also updates t if a step is Chris@16: * performed. Chris@16: * Chris@16: * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Chris@16: * Simple System concept. Chris@16: * \param in The state of the ODE which should be solved. Chris@16: * \param dxdt The derivative of state. Chris@16: * \param t The value of the time. Updated if the step is successful. Chris@16: * \param out Used to store the result of the step. Chris@16: * \param dt The step size. Updated. Chris@16: * \return success if the step was accepted, fail otherwise. Chris@16: */ Chris@16: template< class System , class StateIn , class DerivIn , class StateOut > Chris@16: controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt ) Chris@16: { Chris@16: BOOST_USING_STD_MIN(); Chris@16: BOOST_USING_STD_MAX(); Chris@16: using std::pow; Chris@16: Chris@16: m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ); Chris@16: Chris@16: // do one step with error calculation Chris@16: m_stepper.do_step( system , in , dxdt , t , out , dt , m_xerr.m_v ); Chris@16: Chris@16: m_max_rel_error = m_error_checker.error( m_stepper.algebra() , in , dxdt , m_xerr.m_v , dt ); Chris@16: Chris@16: if( m_max_rel_error > 1.0 ) Chris@16: { Chris@16: // error too large - decrease dt ,limit scaling factor to 0.2 and reset state Chris@101: dt *= max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast( static_cast(9)/static_cast(10) * Chris@101: pow( m_max_rel_error , static_cast(-1) / ( m_stepper.error_order() - 1 ) ) ) , Chris@101: static_cast( static_cast(1)/static_cast (5) ) ); Chris@16: return fail; Chris@16: } Chris@16: else Chris@16: { Chris@16: if( m_max_rel_error < 0.5 ) Chris@16: { Chris@16: // error should be > 0 Chris@101: m_max_rel_error = max BOOST_PREVENT_MACRO_SUBSTITUTION ( Chris@101: static_cast( pow( static_cast(5.0) , -static_cast(m_stepper.stepper_order()) ) ) , Chris@101: m_max_rel_error ); Chris@16: //error too small - increase dt and keep the evolution and limit scaling factor to 5.0 Chris@16: t += dt; Chris@16: dt *= static_cast(9)/static_cast(10) * pow( m_max_rel_error , Chris@16: static_cast(-1) / m_stepper.stepper_order() ); Chris@16: return success; Chris@16: } Chris@16: else Chris@16: { Chris@16: t += dt; Chris@16: return success; Chris@16: } Chris@16: } Chris@16: } Chris@16: Chris@16: /** Chris@16: * \brief Returns the error of the last step. Chris@16: * Chris@16: * returns The last error of the step. Chris@16: */ Chris@16: value_type last_error( void ) const Chris@16: { Chris@16: return m_max_rel_error; Chris@16: } Chris@16: Chris@16: Chris@16: /** Chris@16: * \brief Adjust the size of all temporaries in the stepper manually. Chris@16: * \param x A state from which the size of the temporaries to be resized is deduced. Chris@16: */ Chris@16: template< class StateType > Chris@16: void adjust_size( const StateType &x ) Chris@16: { Chris@16: resize_m_xerr_impl( x ); Chris@16: resize_m_dxdt_impl( x ); Chris@16: resize_m_xnew_impl( x ); Chris@16: m_stepper.adjust_size( x ); Chris@16: } Chris@16: Chris@16: /** Chris@16: * \brief Returns the instance of the underlying stepper. Chris@16: * \returns The instance of the underlying stepper. Chris@16: */ Chris@16: stepper_type& stepper( void ) Chris@16: { Chris@16: return m_stepper; Chris@16: } Chris@16: Chris@16: /** Chris@16: * \brief Returns the instance of the underlying stepper. Chris@16: * \returns The instance of the underlying stepper. Chris@16: */ Chris@16: const stepper_type& stepper( void ) const Chris@16: { Chris@16: return m_stepper; Chris@16: } Chris@16: Chris@16: private: Chris@16: Chris@16: Chris@16: template< class System , class StateInOut > Chris@16: controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt ) Chris@16: { Chris@16: typename odeint::unwrap_reference< System >::type &sys = system; Chris@16: m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); Chris@16: sys( x , m_dxdt.m_v ,t ); Chris@16: return try_step( system , x , m_dxdt.m_v , t , dt ); Chris@16: } Chris@16: Chris@16: template< class StateIn > Chris@16: bool resize_m_xerr_impl( const StateIn &x ) Chris@16: { Chris@16: return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable::type() ); Chris@16: } Chris@16: Chris@16: template< class StateIn > Chris@16: bool resize_m_dxdt_impl( const StateIn &x ) Chris@16: { Chris@16: return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable::type() ); Chris@16: } Chris@16: Chris@16: template< class StateIn > Chris@16: bool resize_m_xnew_impl( const StateIn &x ) Chris@16: { Chris@16: return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable::type() ); Chris@16: } Chris@16: Chris@16: Chris@16: Chris@16: stepper_type m_stepper; Chris@16: error_checker_type m_error_checker; Chris@16: Chris@16: resizer_type m_dxdt_resizer; Chris@16: resizer_type m_xerr_resizer; Chris@16: resizer_type m_xnew_resizer; Chris@16: Chris@16: wrapped_deriv_type m_dxdt; Chris@16: wrapped_state_type m_xerr; Chris@16: wrapped_state_type m_xnew; Chris@16: value_type m_max_rel_error; Chris@16: }; Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: Chris@16: /* Chris@16: * explicit stepper fsal version Chris@16: * Chris@16: * the class introduces the following try_step overloads Chris@16: * try_step( sys , x , t , dt ) Chris@16: * try_step( sys , in , t , out , dt ) Chris@16: * try_step( sys , x , dxdt , t , dt ) Chris@16: * try_step( sys , in , dxdt_in , t , out , dxdt_out , dt ) Chris@16: */ Chris@16: /** Chris@16: * \brief Implements step size control for Runge-Kutta FSAL steppers with Chris@16: * error estimation. Chris@16: * Chris@16: * This class implements the step size control for FSAL Runge-Kutta Chris@16: * steppers with error estimation. Chris@16: * Chris@16: * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept. Chris@16: * \tparam ErrorChecker The error checker Chris@16: * \tparam Resizer The resizer policy type. Chris@16: */ Chris@16: template< Chris@16: class ErrorStepper , Chris@16: class ErrorChecker , Chris@16: class Resizer Chris@16: > Chris@16: class controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_fsal_tag > Chris@16: { Chris@16: Chris@16: public: Chris@16: Chris@16: typedef ErrorStepper stepper_type; Chris@16: typedef typename stepper_type::state_type state_type; Chris@16: typedef typename stepper_type::value_type value_type; Chris@16: typedef typename stepper_type::deriv_type deriv_type; Chris@16: typedef typename stepper_type::time_type time_type; Chris@16: typedef typename stepper_type::algebra_type algebra_type; Chris@16: typedef typename stepper_type::operations_type operations_type; Chris@16: typedef Resizer resizer_type; Chris@16: typedef ErrorChecker error_checker_type; Chris@16: typedef explicit_controlled_stepper_fsal_tag stepper_category; Chris@16: Chris@16: #ifndef DOXYGEN_SKIP Chris@16: typedef typename stepper_type::wrapped_state_type wrapped_state_type; Chris@16: typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type; Chris@16: Chris@16: typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > controlled_stepper_type; Chris@16: #endif // DOXYGEN_SKIP Chris@16: Chris@16: /** Chris@16: * \brief Constructs the controlled Runge-Kutta stepper. Chris@16: * \param error_checker An instance of the error checker. Chris@16: * \param stepper An instance of the underlying stepper. Chris@16: */ Chris@16: controlled_runge_kutta( Chris@16: const error_checker_type &error_checker = error_checker_type() , Chris@16: const stepper_type &stepper = stepper_type() Chris@16: ) Chris@16: : m_stepper( stepper ) , m_error_checker( error_checker ) , Chris@16: m_first_call( true ) Chris@16: { } Chris@16: Chris@16: /* Chris@16: * Version 1 : try_step( sys , x , t , dt ) Chris@16: * Chris@16: * The two overloads are needed in order to solve the forwarding problem Chris@16: */ Chris@16: /** Chris@16: * \brief Tries to perform one step. Chris@16: * Chris@16: * This method tries to do one step with step size dt. If the error estimate Chris@16: * is to large, the step is rejected and the method returns fail and the Chris@16: * step size dt is reduced. If the error estimate is acceptably small, the Chris@16: * step is performed, success is returned and dt might be increased to make Chris@16: * the steps as large as possible. This method also updates t if a step is Chris@16: * performed. Chris@16: * Chris@16: * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Chris@16: * Simple System concept. Chris@16: * \param x The state of the ODE which should be solved. Overwritten if Chris@16: * the step is successful. Chris@16: * \param t The value of the time. Updated if the step is successful. Chris@16: * \param dt The step size. Updated. Chris@16: * \return success if the step was accepted, fail otherwise. Chris@16: */ Chris@16: template< class System , class StateInOut > Chris@16: controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt ) Chris@16: { Chris@16: return try_step_v1( system , x , t , dt ); Chris@16: } Chris@16: Chris@16: Chris@16: /** Chris@16: * \brief Tries to perform one step. Solves the forwarding problem and Chris@16: * allows for using boost range as state_type. Chris@16: * Chris@16: * This method tries to do one step with step size dt. If the error estimate Chris@16: * is to large, the step is rejected and the method returns fail and the Chris@16: * step size dt is reduced. If the error estimate is acceptably small, the Chris@16: * step is performed, success is returned and dt might be increased to make Chris@16: * the steps as large as possible. This method also updates t if a step is Chris@16: * performed. Chris@16: * Chris@16: * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Chris@16: * Simple System concept. Chris@16: * \param x The state of the ODE which should be solved. Overwritten if Chris@16: * the step is successful. Can be a boost range. Chris@16: * \param t The value of the time. Updated if the step is successful. Chris@16: * \param dt The step size. Updated. Chris@16: * \return success if the step was accepted, fail otherwise. Chris@16: */ Chris@16: template< class System , class StateInOut > Chris@16: controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt ) Chris@16: { Chris@16: return try_step_v1( system , x , t , dt ); Chris@16: } Chris@16: Chris@16: Chris@16: Chris@16: /* Chris@16: * Version 2 : try_step( sys , in , t , out , dt ); Chris@16: * Chris@16: * This version does not solve the forwarding problem, boost::range can not be used. Chris@16: * Chris@16: * The disabler is needed to solve ambiguous overloads Chris@16: */ Chris@16: /** Chris@16: * \brief Tries to perform one step. Chris@16: * Chris@16: * \note This method is disabled if state_type=time_type to avoid ambiguity. Chris@16: * Chris@16: * This method tries to do one step with step size dt. If the error estimate Chris@16: * is to large, the step is rejected and the method returns fail and the Chris@16: * step size dt is reduced. If the error estimate is acceptably small, the Chris@16: * step is performed, success is returned and dt might be increased to make Chris@16: * the steps as large as possible. This method also updates t if a step is Chris@16: * performed. Chris@16: * Chris@16: * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Chris@16: * Simple System concept. Chris@16: * \param in The state of the ODE which should be solved. Chris@16: * \param t The value of the time. Updated if the step is successful. Chris@16: * \param out Used to store the result of the step. Chris@16: * \param dt The step size. Updated. Chris@16: * \return success if the step was accepted, fail otherwise. Chris@16: */ Chris@16: template< class System , class StateIn , class StateOut > Chris@16: typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type Chris@16: try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt ) Chris@16: { Chris@16: if( m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) || m_first_call ) Chris@16: { Chris@16: initialize( system , in , t ); Chris@16: } Chris@16: return try_step( system , in , m_dxdt.m_v , t , out , dt ); Chris@16: } Chris@16: Chris@16: Chris@16: /* Chris@16: * Version 3 : try_step( sys , x , dxdt , t , dt ) Chris@16: * Chris@16: * This version does not solve the forwarding problem, boost::range can not be used. Chris@16: */ Chris@16: /** Chris@16: * \brief Tries to perform one step. Chris@16: * Chris@16: * This method tries to do one step with step size dt. If the error estimate Chris@16: * is to large, the step is rejected and the method returns fail and the Chris@16: * step size dt is reduced. If the error estimate is acceptably small, the Chris@16: * step is performed, success is returned and dt might be increased to make Chris@16: * the steps as large as possible. This method also updates t if a step is Chris@16: * performed. Chris@16: * Chris@16: * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Chris@16: * Simple System concept. Chris@16: * \param x The state of the ODE which should be solved. Overwritten if Chris@16: * the step is successful. Chris@16: * \param dxdt The derivative of state. Chris@16: * \param t The value of the time. Updated if the step is successful. Chris@16: * \param dt The step size. Updated. Chris@16: * \return success if the step was accepted, fail otherwise. Chris@16: */ Chris@16: template< class System , class StateInOut , class DerivInOut > Chris@16: controlled_step_result try_step( System system , StateInOut &x , DerivInOut &dxdt , time_type &t , time_type &dt ) Chris@16: { Chris@16: m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); Chris@16: m_dxdt_new_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_new_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); Chris@16: controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , m_dxdtnew.m_v , dt ); Chris@16: if( res == success ) Chris@16: { Chris@16: boost::numeric::odeint::copy( m_xnew.m_v , x ); Chris@16: boost::numeric::odeint::copy( m_dxdtnew.m_v , dxdt ); Chris@16: } Chris@16: return res; Chris@16: } Chris@16: Chris@16: Chris@16: /* Chris@16: * Version 4 : try_step( sys , in , dxdt_in , t , out , dxdt_out , dt ) Chris@16: * Chris@16: * This version does not solve the forwarding problem, boost::range can not be used. Chris@16: */ Chris@16: /** Chris@16: * \brief Tries to perform one step. Chris@16: * Chris@16: * This method tries to do one step with step size dt. If the error estimate Chris@16: * is to large, the step is rejected and the method returns fail and the Chris@16: * step size dt is reduced. If the error estimate is acceptably small, the Chris@16: * step is performed, success is returned and dt might be increased to make Chris@16: * the steps as large as possible. This method also updates t if a step is Chris@16: * performed. Chris@16: * Chris@16: * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Chris@16: * Simple System concept. Chris@16: * \param in The state of the ODE which should be solved. Chris@16: * \param dxdt The derivative of state. Chris@16: * \param t The value of the time. Updated if the step is successful. Chris@16: * \param out Used to store the result of the step. Chris@16: * \param dt The step size. Updated. Chris@16: * \return success if the step was accepted, fail otherwise. Chris@16: */ Chris@16: template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut > Chris@16: controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t , Chris@16: StateOut &out , DerivOut &dxdt_out , time_type &dt ) Chris@16: { Chris@16: BOOST_USING_STD_MIN(); Chris@16: BOOST_USING_STD_MAX(); Chris@16: Chris@16: using std::pow; Chris@16: Chris@16: m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ); Chris@16: Chris@16: //fsal: m_stepper.get_dxdt( dxdt ); Chris@16: //fsal: m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err ); Chris@16: m_stepper.do_step( system , in , dxdt_in , t , out , dxdt_out , dt , m_xerr.m_v ); Chris@16: Chris@16: // this potentially overwrites m_x_err! (standard_error_checker does, at least) Chris@16: value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt_in , m_xerr.m_v , dt ); Chris@16: Chris@16: if( max_rel_err > 1.0 ) Chris@16: { Chris@16: // error too large - decrease dt ,limit scaling factor to 0.2 and reset state Chris@101: dt *= max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast( static_cast(9)/static_cast(10) * Chris@101: pow( max_rel_err , static_cast(-1) / ( m_stepper.error_order() - 1 ) ) ) , Chris@101: static_cast( static_cast(1)/static_cast (5)) ); Chris@16: return fail; Chris@16: } Chris@16: else Chris@16: { Chris@16: if( max_rel_err < 0.5 ) Chris@16: { //error too small - increase dt and keep the evolution and limit scaling factor to 5.0 Chris@16: // error should be > 0 Chris@101: max_rel_err = max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast( pow( static_cast(5.0) , -static_cast(m_stepper.stepper_order()) ) ) , Chris@101: max_rel_err ); Chris@16: t += dt; Chris@16: dt *= static_cast( static_cast(9)/static_cast(10) * pow( max_rel_err , static_cast(-1) / m_stepper.stepper_order() ) ); Chris@16: return success; Chris@16: } Chris@16: else Chris@16: { Chris@16: t += dt; Chris@16: return success; Chris@16: } Chris@16: } Chris@16: } Chris@16: Chris@16: Chris@16: /** Chris@16: * \brief Resets the internal state of the underlying FSAL stepper. Chris@16: */ Chris@16: void reset( void ) Chris@16: { Chris@16: m_first_call = true; Chris@16: } Chris@16: Chris@16: /** Chris@16: * \brief Initializes the internal state storing an internal copy of the derivative. Chris@16: * Chris@16: * \param deriv The initial derivative of the ODE. Chris@16: */ Chris@16: template< class DerivIn > Chris@16: void initialize( const DerivIn &deriv ) Chris@16: { Chris@16: boost::numeric::odeint::copy( deriv , m_dxdt.m_v ); Chris@16: m_first_call = false; Chris@16: } Chris@16: Chris@16: /** Chris@16: * \brief Initializes the internal state storing an internal copy of the derivative. Chris@16: * Chris@16: * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Chris@16: * Simple System concept. Chris@16: * \param x The initial state of the ODE which should be solved. Chris@16: * \param t The initial time. Chris@16: */ Chris@16: template< class System , class StateIn > Chris@16: void initialize( System system , const StateIn &x , time_type t ) Chris@16: { Chris@16: typename odeint::unwrap_reference< System >::type &sys = system; Chris@16: sys( x , m_dxdt.m_v , t ); Chris@16: m_first_call = false; Chris@16: } Chris@16: Chris@16: /** Chris@16: * \brief Returns true if the stepper has been initialized, false otherwise. Chris@16: * Chris@16: * \return true, if the stepper has been initialized, false otherwise. Chris@16: */ Chris@16: bool is_initialized( void ) const Chris@16: { Chris@16: return ! m_first_call; Chris@16: } Chris@16: Chris@16: Chris@16: /** Chris@16: * \brief Adjust the size of all temporaries in the stepper manually. Chris@16: * \param x A state from which the size of the temporaries to be resized is deduced. Chris@16: */ Chris@16: template< class StateType > Chris@16: void adjust_size( const StateType &x ) Chris@16: { Chris@16: resize_m_xerr_impl( x ); Chris@16: resize_m_dxdt_impl( x ); Chris@16: resize_m_dxdt_new_impl( x ); Chris@16: resize_m_xnew_impl( x ); Chris@16: } Chris@16: Chris@16: Chris@16: /** Chris@16: * \brief Returns the instance of the underlying stepper. Chris@16: * \returns The instance of the underlying stepper. Chris@16: */ Chris@16: stepper_type& stepper( void ) Chris@16: { Chris@16: return m_stepper; Chris@16: } Chris@16: Chris@16: /** Chris@16: * \brief Returns the instance of the underlying stepper. Chris@16: * \returns The instance of the underlying stepper. Chris@16: */ Chris@16: const stepper_type& stepper( void ) const Chris@16: { Chris@16: return m_stepper; Chris@16: } Chris@16: Chris@16: Chris@16: Chris@16: private: Chris@16: Chris@16: Chris@16: template< class StateIn > Chris@16: bool resize_m_xerr_impl( const StateIn &x ) Chris@16: { Chris@16: return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable::type() ); Chris@16: } Chris@16: Chris@16: template< class StateIn > Chris@16: bool resize_m_dxdt_impl( const StateIn &x ) Chris@16: { Chris@16: return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable::type() ); Chris@16: } Chris@16: Chris@16: template< class StateIn > Chris@16: bool resize_m_dxdt_new_impl( const StateIn &x ) Chris@16: { Chris@16: return adjust_size_by_resizeability( m_dxdtnew , x , typename is_resizeable::type() ); Chris@16: } Chris@16: Chris@16: template< class StateIn > Chris@16: bool resize_m_xnew_impl( const StateIn &x ) Chris@16: { Chris@16: return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable::type() ); Chris@16: } Chris@16: Chris@16: Chris@16: template< class System , class StateInOut > Chris@16: controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt ) Chris@16: { Chris@16: if( m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ) || m_first_call ) Chris@16: { Chris@16: initialize( system , x , t ); Chris@16: } Chris@16: return try_step( system , x , m_dxdt.m_v , t , dt ); Chris@16: } Chris@16: Chris@16: Chris@16: stepper_type m_stepper; Chris@16: error_checker_type m_error_checker; Chris@16: Chris@16: resizer_type m_dxdt_resizer; Chris@16: resizer_type m_xerr_resizer; Chris@16: resizer_type m_xnew_resizer; Chris@16: resizer_type m_dxdt_new_resizer; Chris@16: Chris@16: wrapped_deriv_type m_dxdt; Chris@16: wrapped_state_type m_xerr; Chris@16: wrapped_state_type m_xnew; Chris@16: wrapped_deriv_type m_dxdtnew; Chris@16: bool m_first_call; Chris@16: }; Chris@16: Chris@16: Chris@16: /********** DOXYGEN **********/ Chris@16: Chris@16: /**** DEFAULT ERROR CHECKER ****/ Chris@16: Chris@16: /** Chris@16: * \class default_error_checker Chris@16: * \brief The default error checker to be used with Runge-Kutta error steppers Chris@16: * Chris@16: * This class provides the default mechanism to compare the error estimates Chris@16: * reported by Runge-Kutta error steppers with user defined error bounds. Chris@16: * It is used by the controlled_runge_kutta steppers. Chris@16: * Chris@16: * \tparam Value The value type. Chris@16: * \tparam Algebra The algebra type. Chris@16: * \tparam Operations The operations type. Chris@16: */ Chris@16: Chris@16: /** Chris@16: * \fn default_error_checker( value_type eps_abs , value_type eps_rel , value_type a_x , value_type a_dxdt ) Chris@16: * \brief Constructs the error checker. Chris@16: * Chris@16: * The error is calculated as follows: ???? Chris@16: * Chris@16: * \param eps_abs Absolute tolerance level. Chris@16: * \param eps_rel Relative tolerance level. Chris@16: * \param a_x Factor for the weight of the state. Chris@16: * \param a_dxdt Factor for the weight of the derivative. Chris@16: */ Chris@16: Chris@16: /** Chris@16: * \fn error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const Chris@16: * \brief Calculates the error level. Chris@16: * Chris@16: * If the returned error level is greater than 1, the estimated error was Chris@16: * larger than the permitted error bounds and the step should be repeated Chris@16: * with a smaller step size. Chris@16: * Chris@16: * \param x_old State at the beginning of the step. Chris@16: * \param dxdt_old Derivative at the beginning of the step. Chris@16: * \param x_err Error estimate. Chris@16: * \param dt Time step. Chris@16: * \return error Chris@16: */ Chris@16: Chris@16: /** Chris@16: * \fn error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const Chris@16: * \brief Calculates the error level using a given algebra. Chris@16: * Chris@16: * If the returned error level is greater than 1, the estimated error was Chris@16: * larger than the permitted error bounds and the step should be repeated Chris@16: * with a smaller step size. Chris@16: * Chris@16: * \param algebra The algebra used for calculation of the error. Chris@16: * \param x_old State at the beginning of the step. Chris@16: * \param dxdt_old Derivative at the beginning of the step. Chris@16: * \param x_err Error estimate. Chris@16: * \param dt Time step. Chris@16: * \return error Chris@16: */ Chris@16: Chris@16: Chris@16: } // odeint Chris@16: } // numeric Chris@16: } // boost Chris@16: Chris@16: Chris@16: #endif // BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED