| /* [auto_generated] |
| boost/numeric/odeint/stepper/controlled_runge_kutta.hpp |
| |
| [begin_description] |
| The default controlled stepper which can be used with all explicit Runge-Kutta error steppers. |
| [end_description] |
| |
| Copyright 2010-2013 Karsten Ahnert |
| Copyright 2010-2013 Mario Mulansky |
| Copyright 2012 Christoph Koke |
| |
| 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_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED |
| #define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED |
| |
| |
| |
| #include <cmath> |
| |
| #include <boost/config.hpp> |
| #include <boost/utility/enable_if.hpp> |
| #include <boost/type_traits/is_same.hpp> |
| |
| #include <boost/numeric/odeint/util/bind.hpp> |
| #include <boost/numeric/odeint/util/unwrap_reference.hpp> |
| #include <boost/numeric/odeint/util/copy.hpp> |
| |
| #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/algebra/range_algebra.hpp> |
| #include <boost/numeric/odeint/algebra/default_operations.hpp> |
| #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> |
| |
| #include <boost/numeric/odeint/stepper/controlled_step_result.hpp> |
| #include <boost/numeric/odeint/stepper/stepper_categories.hpp> |
| |
| namespace boost { |
| namespace numeric { |
| namespace odeint { |
| |
| |
| template |
| < |
| class Value , |
| class Algebra , |
| class Operations |
| > |
| class default_error_checker |
| { |
| public: |
| |
| typedef Value value_type; |
| typedef Algebra algebra_type; |
| typedef Operations operations_type; |
| |
| default_error_checker( |
| 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 ) |
| { } |
| |
| |
| 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 > |
| value_type error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const |
| { |
| // 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 ) ) ); |
| |
| // value_type res = algebra.reduce( x_err , |
| // typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0 ) ); |
| return algebra.norm_inf( x_err ); |
| } |
| |
| private: |
| |
| value_type m_eps_abs; |
| value_type m_eps_rel; |
| value_type m_a_x; |
| value_type m_a_dxdt; |
| |
| }; |
| |
| |
| |
| |
| |
| |
| |
| |
| /* |
| * error stepper category dispatcher |
| */ |
| template< |
| class ErrorStepper , |
| class ErrorChecker = default_error_checker< typename ErrorStepper::value_type , |
| typename ErrorStepper::algebra_type , |
| typename ErrorStepper::operations_type > , |
| class Resizer = typename ErrorStepper::resizer_type , |
| class ErrorStepperCategory = typename ErrorStepper::stepper_category |
| > |
| class controlled_runge_kutta ; |
| |
| |
| |
| /* |
| * explicit stepper version |
| * |
| * this class introduces the following try_step overloads |
| * try_step( sys , x , t , dt ) |
| * try_step( sys , x , dxdt , t , dt ) |
| * try_step( sys , in , t , out , dt ) |
| * try_step( sys , in , dxdt , t , out , dt ) |
| */ |
| /** |
| * \brief Implements step size control for Runge-Kutta steppers with error |
| * estimation. |
| * |
| * This class implements the step size control for standard Runge-Kutta |
| * steppers with error estimation. |
| * |
| * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept. |
| * \tparam ErrorChecker The error checker |
| * \tparam Resizer The resizer policy type. |
| */ |
| template< |
| class ErrorStepper , |
| class ErrorChecker , |
| class Resizer |
| > |
| class controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > |
| { |
| |
| public: |
| |
| typedef ErrorStepper stepper_type; |
| 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; |
| typedef typename stepper_type::algebra_type algebra_type; |
| typedef typename stepper_type::operations_type operations_type; |
| typedef Resizer resizer_type; |
| typedef ErrorChecker error_checker_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; |
| #endif //DOXYGEN_SKIP |
| |
| |
| /** |
| * \brief Constructs the controlled Runge-Kutta stepper. |
| * \param error_checker An instance of the error checker. |
| * \param stepper An instance of the underlying stepper. |
| */ |
| controlled_runge_kutta( |
| const error_checker_type &error_checker = error_checker_type( ) , |
| const stepper_type &stepper = stepper_type( ) |
| ) |
| : m_stepper( stepper ) , m_error_checker( error_checker ) |
| { } |
| |
| |
| |
| /* |
| * Version 1 : try_step( sys , x , t , dt ) |
| * |
| * The overloads are needed to solve the forwarding problem |
| */ |
| /** |
| * \brief Tries to perform one step. |
| * |
| * This method tries to do one step with step size dt. If the error estimate |
| * is to large, the step is rejected and the method returns fail and the |
| * step size dt is reduced. If the error estimate is acceptably small, the |
| * step is performed, success is returned and dt might be increased to make |
| * the steps as large as possible. This method also updates t if a step is |
| * performed. |
| * |
| * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the |
| * Simple System concept. |
| * \param x The state of the ODE which should be solved. Overwritten if |
| * the step is successful. |
| * \param t The value of the time. Updated if the step is successful. |
| * \param dt The step size. Updated. |
| * \return success if the step was accepted, fail otherwise. |
| */ |
| template< class System , class StateInOut > |
| controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt ) |
| { |
| return try_step_v1( system , x , t, dt ); |
| } |
| |
| /** |
| * \brief Tries to perform one step. Solves the forwarding problem and |
| * allows for using boost range as state_type. |
| * |
| * This method tries to do one step with step size dt. If the error estimate |
| * is to large, the step is rejected and the method returns fail and the |
| * step size dt is reduced. If the error estimate is acceptably small, the |
| * step is performed, success is returned and dt might be increased to make |
| * the steps as large as possible. This method also updates t if a step is |
| * performed. |
| * |
| * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the |
| * Simple System concept. |
| * \param x The state of the ODE which should be solved. Overwritten if |
| * the step is successful. Can be a boost range. |
| * \param t The value of the time. Updated if the step is successful. |
| * \param dt The step size. Updated. |
| * \return success if the step was accepted, fail otherwise. |
| */ |
| template< class System , class StateInOut > |
| controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt ) |
| { |
| return try_step_v1( system , x , t, dt ); |
| } |
| |
| |
| |
| /* |
| * Version 2 : try_step( sys , x , dxdt , t , dt ) |
| * |
| * this version does not solve the forwarding problem, boost.range can not be used |
| */ |
| /** |
| * \brief Tries to perform one step. |
| * |
| * This method tries to do one step with step size dt. If the error estimate |
| * is to large, the step is rejected and the method returns fail and the |
| * step size dt is reduced. If the error estimate is acceptably small, the |
| * step is performed, success is returned and dt might be increased to make |
| * the steps as large as possible. This method also updates t if a step is |
| * performed. |
| * |
| * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the |
| * Simple System concept. |
| * \param x The state of the ODE which should be solved. Overwritten if |
| * the step is successful. |
| * \param dxdt The derivative of state. |
| * \param t The value of the time. Updated if the step is successful. |
| * \param dt The step size. Updated. |
| * \return success if the step was accepted, fail otherwise. |
| */ |
| template< class System , class StateInOut , class DerivIn > |
| controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt ) |
| { |
| m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); |
| controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt ); |
| if( res == success ) |
| { |
| boost::numeric::odeint::copy( m_xnew.m_v , x ); |
| } |
| return res; |
| } |
| |
| /* |
| * Version 3 : try_step( sys , in , t , out , dt ) |
| * |
| * this version does not solve the forwarding problem, boost.range can not be used |
| * |
| * the disable is needed to avoid ambiguous overloads if state_type = time_type |
| */ |
| /** |
| * \brief Tries to perform one step. |
| * |
| * \note This method is disabled if state_type=time_type to avoid ambiguity. |
| * |
| * This method tries to do one step with step size dt. If the error estimate |
| * is to large, the step is rejected and the method returns fail and the |
| * step size dt is reduced. If the error estimate is acceptably small, the |
| * step is performed, success is returned and dt might be increased to make |
| * the steps as large as possible. This method also updates t if a step is |
| * performed. |
| * |
| * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the |
| * Simple System concept. |
| * \param in The state of the ODE which should be solved. |
| * \param t The value of the time. Updated if the step is successful. |
| * \param out Used to store the result of the step. |
| * \param dt The step size. Updated. |
| * \return success if the step was accepted, fail otherwise. |
| */ |
| template< class System , class StateIn , class StateOut > |
| typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type |
| try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt ) |
| { |
| typename odeint::unwrap_reference< System >::type &sys = system; |
| m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ); |
| sys( in , m_dxdt.m_v , t ); |
| return try_step( system , in , m_dxdt.m_v , t , out , dt ); |
| } |
| |
| |
| /* |
| * Version 4 : try_step( sys , in , dxdt , t , out , dt ) |
| * |
| * this version does not solve the forwarding problem, boost.range can not be used |
| */ |
| /** |
| * \brief Tries to perform one step. |
| * |
| * This method tries to do one step with step size dt. If the error estimate |
| * is to large, the step is rejected and the method returns fail and the |
| * step size dt is reduced. If the error estimate is acceptably small, the |
| * step is performed, success is returned and dt might be increased to make |
| * the steps as large as possible. This method also updates t if a step is |
| * performed. |
| * |
| * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the |
| * Simple System concept. |
| * \param in The state of the ODE which should be solved. |
| * \param dxdt The derivative of state. |
| * \param t The value of the time. Updated if the step is successful. |
| * \param out Used to store the result of the step. |
| * \param dt The step size. Updated. |
| * \return success if the step was accepted, fail otherwise. |
| */ |
| 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; |
| |
| 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 ); |
| |
| if( m_max_rel_error > 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) ) ); |
| return fail; |
| } |
| 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; |
| } |
| } |
| } |
| |
| /** |
| * \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. |
| */ |
| template< class StateType > |
| void adjust_size( const StateType &x ) |
| { |
| resize_m_xerr_impl( x ); |
| resize_m_dxdt_impl( x ); |
| resize_m_xnew_impl( x ); |
| m_stepper.adjust_size( x ); |
| } |
| |
| /** |
| * \brief Returns the instance of the underlying stepper. |
| * \returns The instance of the underlying stepper. |
| */ |
| stepper_type& stepper( void ) |
| { |
| return m_stepper; |
| } |
| |
| /** |
| * \brief Returns the instance of the underlying stepper. |
| * \returns The instance of the underlying stepper. |
| */ |
| const stepper_type& stepper( void ) const |
| { |
| return m_stepper; |
| } |
| |
| private: |
| |
| |
| template< class System , class StateInOut > |
| controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt ) |
| { |
| typename odeint::unwrap_reference< System >::type &sys = system; |
| m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); |
| sys( x , m_dxdt.m_v ,t ); |
| return try_step( system , x , m_dxdt.m_v , t , dt ); |
| } |
| |
| template< class StateIn > |
| bool resize_m_xerr_impl( const StateIn &x ) |
| { |
| return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() ); |
| } |
| |
| template< class StateIn > |
| bool resize_m_dxdt_impl( const StateIn &x ) |
| { |
| return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() ); |
| } |
| |
| template< class StateIn > |
| bool resize_m_xnew_impl( const StateIn &x ) |
| { |
| return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() ); |
| } |
| |
| |
| |
| stepper_type m_stepper; |
| error_checker_type m_error_checker; |
| |
| resizer_type m_dxdt_resizer; |
| resizer_type m_xerr_resizer; |
| resizer_type m_xnew_resizer; |
| |
| wrapped_deriv_type m_dxdt; |
| wrapped_state_type m_xerr; |
| wrapped_state_type m_xnew; |
| value_type m_max_rel_error; |
| }; |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| /* |
| * explicit stepper fsal version |
| * |
| * the class introduces the following try_step overloads |
| * try_step( sys , x , t , dt ) |
| * try_step( sys , in , t , out , dt ) |
| * try_step( sys , x , dxdt , t , dt ) |
| * try_step( sys , in , dxdt_in , t , out , dxdt_out , dt ) |
| */ |
| /** |
| * \brief Implements step size control for Runge-Kutta FSAL steppers with |
| * error estimation. |
| * |
| * This class implements the step size control for FSAL Runge-Kutta |
| * steppers with error estimation. |
| * |
| * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept. |
| * \tparam ErrorChecker The error checker |
| * \tparam Resizer The resizer policy type. |
| */ |
| template< |
| class ErrorStepper , |
| class ErrorChecker , |
| class Resizer |
| > |
| class controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_fsal_tag > |
| { |
| |
| public: |
| |
| typedef ErrorStepper stepper_type; |
| 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; |
| typedef typename stepper_type::algebra_type algebra_type; |
| typedef typename stepper_type::operations_type operations_type; |
| typedef Resizer resizer_type; |
| typedef ErrorChecker error_checker_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; |
| #endif // DOXYGEN_SKIP |
| |
| /** |
| * \brief Constructs the controlled Runge-Kutta stepper. |
| * \param error_checker An instance of the error checker. |
| * \param stepper An instance of the underlying stepper. |
| */ |
| controlled_runge_kutta( |
| const error_checker_type &error_checker = error_checker_type() , |
| const stepper_type &stepper = stepper_type() |
| ) |
| : m_stepper( stepper ) , m_error_checker( error_checker ) , |
| m_first_call( true ) |
| { } |
| |
| /* |
| * Version 1 : try_step( sys , x , t , dt ) |
| * |
| * The two overloads are needed in order to solve the forwarding problem |
| */ |
| /** |
| * \brief Tries to perform one step. |
| * |
| * This method tries to do one step with step size dt. If the error estimate |
| * is to large, the step is rejected and the method returns fail and the |
| * step size dt is reduced. If the error estimate is acceptably small, the |
| * step is performed, success is returned and dt might be increased to make |
| * the steps as large as possible. This method also updates t if a step is |
| * performed. |
| * |
| * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the |
| * Simple System concept. |
| * \param x The state of the ODE which should be solved. Overwritten if |
| * the step is successful. |
| * \param t The value of the time. Updated if the step is successful. |
| * \param dt The step size. Updated. |
| * \return success if the step was accepted, fail otherwise. |
| */ |
| template< class System , class StateInOut > |
| controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt ) |
| { |
| return try_step_v1( system , x , t , dt ); |
| } |
| |
| |
| /** |
| * \brief Tries to perform one step. Solves the forwarding problem and |
| * allows for using boost range as state_type. |
| * |
| * This method tries to do one step with step size dt. If the error estimate |
| * is to large, the step is rejected and the method returns fail and the |
| * step size dt is reduced. If the error estimate is acceptably small, the |
| * step is performed, success is returned and dt might be increased to make |
| * the steps as large as possible. This method also updates t if a step is |
| * performed. |
| * |
| * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the |
| * Simple System concept. |
| * \param x The state of the ODE which should be solved. Overwritten if |
| * the step is successful. Can be a boost range. |
| * \param t The value of the time. Updated if the step is successful. |
| * \param dt The step size. Updated. |
| * \return success if the step was accepted, fail otherwise. |
| */ |
| template< class System , class StateInOut > |
| controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt ) |
| { |
| return try_step_v1( system , x , t , dt ); |
| } |
| |
| |
| |
| /* |
| * Version 2 : try_step( sys , in , t , out , dt ); |
| * |
| * This version does not solve the forwarding problem, boost::range can not be used. |
| * |
| * The disabler is needed to solve ambiguous overloads |
| */ |
| /** |
| * \brief Tries to perform one step. |
| * |
| * \note This method is disabled if state_type=time_type to avoid ambiguity. |
| * |
| * This method tries to do one step with step size dt. If the error estimate |
| * is to large, the step is rejected and the method returns fail and the |
| * step size dt is reduced. If the error estimate is acceptably small, the |
| * step is performed, success is returned and dt might be increased to make |
| * the steps as large as possible. This method also updates t if a step is |
| * performed. |
| * |
| * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the |
| * Simple System concept. |
| * \param in The state of the ODE which should be solved. |
| * \param t The value of the time. Updated if the step is successful. |
| * \param out Used to store the result of the step. |
| * \param dt The step size. Updated. |
| * \return success if the step was accepted, fail otherwise. |
| */ |
| template< class System , class StateIn , class StateOut > |
| typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type |
| try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt ) |
| { |
| 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 ) |
| { |
| initialize( system , in , t ); |
| } |
| return try_step( system , in , m_dxdt.m_v , t , out , dt ); |
| } |
| |
| |
| /* |
| * Version 3 : try_step( sys , x , dxdt , t , dt ) |
| * |
| * This version does not solve the forwarding problem, boost::range can not be used. |
| */ |
| /** |
| * \brief Tries to perform one step. |
| * |
| * This method tries to do one step with step size dt. If the error estimate |
| * is to large, the step is rejected and the method returns fail and the |
| * step size dt is reduced. If the error estimate is acceptably small, the |
| * step is performed, success is returned and dt might be increased to make |
| * the steps as large as possible. This method also updates t if a step is |
| * performed. |
| * |
| * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the |
| * Simple System concept. |
| * \param x The state of the ODE which should be solved. Overwritten if |
| * the step is successful. |
| * \param dxdt The derivative of state. |
| * \param t The value of the time. Updated if the step is successful. |
| * \param dt The step size. Updated. |
| * \return success if the step was accepted, fail otherwise. |
| */ |
| template< class System , class StateInOut , class DerivInOut > |
| controlled_step_result try_step( System system , StateInOut &x , DerivInOut &dxdt , time_type &t , time_type &dt ) |
| { |
| m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); |
| m_dxdt_new_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_new_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); |
| controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , m_dxdtnew.m_v , dt ); |
| if( res == success ) |
| { |
| boost::numeric::odeint::copy( m_xnew.m_v , x ); |
| boost::numeric::odeint::copy( m_dxdtnew.m_v , dxdt ); |
| } |
| return res; |
| } |
| |
| |
| /* |
| * Version 4 : try_step( sys , in , dxdt_in , t , out , dxdt_out , dt ) |
| * |
| * This version does not solve the forwarding problem, boost::range can not be used. |
| */ |
| /** |
| * \brief Tries to perform one step. |
| * |
| * This method tries to do one step with step size dt. If the error estimate |
| * is to large, the step is rejected and the method returns fail and the |
| * step size dt is reduced. If the error estimate is acceptably small, the |
| * step is performed, success is returned and dt might be increased to make |
| * the steps as large as possible. This method also updates t if a step is |
| * performed. |
| * |
| * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the |
| * Simple System concept. |
| * \param in The state of the ODE which should be solved. |
| * \param dxdt The derivative of state. |
| * \param t The value of the time. Updated if the step is successful. |
| * \param out Used to store the result of the step. |
| * \param dt The step size. Updated. |
| * \return success if the step was accepted, fail otherwise. |
| */ |
| template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut > |
| 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; |
| |
| m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ); |
| |
| //fsal: m_stepper.get_dxdt( dxdt ); |
| //fsal: m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err ); |
| m_stepper.do_step( system , in , dxdt_in , t , out , dxdt_out , dt , m_xerr.m_v ); |
| |
| // this potentially overwrites m_x_err! (standard_error_checker does, at least) |
| value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt_in , m_xerr.m_v , dt ); |
| |
| 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)) ); |
| 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; |
| } |
| } |
| } |
| |
| |
| /** |
| * \brief Resets the internal state of the underlying FSAL stepper. |
| */ |
| void reset( void ) |
| { |
| m_first_call = true; |
| } |
| |
| /** |
| * \brief Initializes the internal state storing an internal copy of the derivative. |
| * |
| * \param deriv The initial derivative of the ODE. |
| */ |
| template< class DerivIn > |
| void initialize( const DerivIn &deriv ) |
| { |
| boost::numeric::odeint::copy( deriv , m_dxdt.m_v ); |
| m_first_call = false; |
| } |
| |
| /** |
| * \brief Initializes the internal state storing an internal copy of the derivative. |
| * |
| * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the |
| * Simple System concept. |
| * \param x The initial state of the ODE which should be solved. |
| * \param t The initial time. |
| */ |
| template< class System , class StateIn > |
| void initialize( System system , const StateIn &x , time_type t ) |
| { |
| typename odeint::unwrap_reference< System >::type &sys = system; |
| sys( x , m_dxdt.m_v , t ); |
| m_first_call = false; |
| } |
| |
| /** |
| * \brief Returns true if the stepper has been initialized, false otherwise. |
| * |
| * \return true, if the stepper has been initialized, false otherwise. |
| */ |
| bool is_initialized( void ) const |
| { |
| return ! m_first_call; |
| } |
| |
| |
| /** |
| * \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. |
| */ |
| template< class StateType > |
| void adjust_size( const StateType &x ) |
| { |
| resize_m_xerr_impl( x ); |
| resize_m_dxdt_impl( x ); |
| resize_m_dxdt_new_impl( x ); |
| resize_m_xnew_impl( x ); |
| } |
| |
| |
| /** |
| * \brief Returns the instance of the underlying stepper. |
| * \returns The instance of the underlying stepper. |
| */ |
| stepper_type& stepper( void ) |
| { |
| return m_stepper; |
| } |
| |
| /** |
| * \brief Returns the instance of the underlying stepper. |
| * \returns The instance of the underlying stepper. |
| */ |
| const stepper_type& stepper( void ) const |
| { |
| return m_stepper; |
| } |
| |
| |
| |
| private: |
| |
| |
| template< class StateIn > |
| bool resize_m_xerr_impl( const StateIn &x ) |
| { |
| return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() ); |
| } |
| |
| template< class StateIn > |
| bool resize_m_dxdt_impl( const StateIn &x ) |
| { |
| return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() ); |
| } |
| |
| template< class StateIn > |
| bool resize_m_dxdt_new_impl( const StateIn &x ) |
| { |
| return adjust_size_by_resizeability( m_dxdtnew , x , typename is_resizeable<deriv_type>::type() ); |
| } |
| |
| template< class StateIn > |
| bool resize_m_xnew_impl( const StateIn &x ) |
| { |
| return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() ); |
| } |
| |
| |
| template< class System , class StateInOut > |
| controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt ) |
| { |
| 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 ) |
| { |
| initialize( system , x , t ); |
| } |
| return try_step( system , x , m_dxdt.m_v , t , dt ); |
| } |
| |
| |
| stepper_type m_stepper; |
| error_checker_type m_error_checker; |
| |
| resizer_type m_dxdt_resizer; |
| resizer_type m_xerr_resizer; |
| resizer_type m_xnew_resizer; |
| resizer_type m_dxdt_new_resizer; |
| |
| wrapped_deriv_type m_dxdt; |
| wrapped_state_type m_xerr; |
| wrapped_state_type m_xnew; |
| wrapped_deriv_type m_dxdtnew; |
| bool m_first_call; |
| }; |
| |
| |
| /********** DOXYGEN **********/ |
| |
| /**** DEFAULT ERROR CHECKER ****/ |
| |
| /** |
| * \class default_error_checker |
| * \brief The default error checker to be used with Runge-Kutta error steppers |
| * |
| * This class provides the default mechanism to compare the error estimates |
| * reported by Runge-Kutta error steppers with user defined error bounds. |
| * It is used by the controlled_runge_kutta steppers. |
| * |
| * \tparam Value The value 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 ) |
| * \brief Constructs the error checker. |
| * |
| * The error is calculated as follows: ???? |
| * |
| * \param eps_abs Absolute tolerance level. |
| * \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. |
| */ |
| |
| /** |
| * \fn error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const |
| * \brief Calculates the error level. |
| * |
| * If the returned error level is greater than 1, the estimated error was |
| * larger than the permitted error bounds and the step should be repeated |
| * with a smaller step size. |
| * |
| * \param x_old State at the beginning of the step. |
| * \param dxdt_old Derivative at the beginning of the step. |
| * \param x_err Error estimate. |
| * \param dt Time step. |
| * \return error |
| */ |
| |
| /** |
| * \fn error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const |
| * \brief Calculates the error level using a given algebra. |
| * |
| * If the returned error level is greater than 1, the estimated error was |
| * larger than the permitted error bounds and the step should be repeated |
| * with a smaller step size. |
| * |
| * \param algebra The algebra used for calculation of the error. |
| * \param x_old State at the beginning of the step. |
| * \param dxdt_old Derivative at the beginning of the step. |
| * \param x_err Error estimate. |
| * \param dt Time step. |
| * \return error |
| */ |
| |
| |
| } // odeint |
| } // numeric |
| } // boost |
| |
| |
| #endif // BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED |
