odeint/stepper/controlled_runge_kutta.hpp

0001 /* [auto_generated]
0002  boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
0003
0004  [begin_description]
0005  The default controlled stepper which can be used with all explicit Runge-Kutta error steppers.
0006  [end_description]
0007
0008  Copyright 2010-2013 Karsten Ahnert
0009  Copyright 2010-2015 Mario Mulansky
0010  Copyright 2012 Christoph Koke
0011
0012  Distributed under the Boost Software License, Version 1.0.
0013  (See accompanying file LICENSE_1_0.txt or
0014  copy at http://www.boost.org/LICENSE_1_0.txt)
0015  */
0016
0017
0018 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED
0019 #define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED
0020
0021
0022
0023 #include <cmath>
0024
0025 #include <boost/config.hpp>
0026 #include <boost/utility/enable_if.hpp>
0027 #include <boost/type_traits/is_same.hpp>
0028
0029 #include <boost/numeric/odeint/util/bind.hpp>
0030 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0031 #include <boost/numeric/odeint/util/copy.hpp>
0032
0033 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0034 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0035 #include <boost/numeric/odeint/util/resizer.hpp>
0036 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
0037
0038 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
0039 #include <boost/numeric/odeint/algebra/default_operations.hpp>
0040 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
0041
0042 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
0043 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0044
0045 namespace boost {
0046 namespace numeric {
0047 namespace odeint {
0048
0049
0050 template
0051 <
0052 class Value ,
0053 class Algebra ,
0054 class Operations
0055 >
0056 class default_error_checker
0057 {
0058 public:
0059
0060     typedef Value value_type;
0061     typedef Algebra algebra_type;
0062     typedef Operations operations_type;
0063
0064     default_error_checker(
0065             value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
0066             value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
0067             value_type a_x = static_cast< value_type >( 1 ) ,
0068             value_type a_dxdt = static_cast< value_type >( 1 ))
0069         : m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
0070     { }
0071
0072
0073     template< class State , class Deriv , class Err, class Time >
0074     value_type error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
0075     {
0076         return error( algebra_type() , x_old , dxdt_old , x_err , dt );
0077     }
0078
0079     template< class State , class Deriv , class Err, class Time >
0080     value_type error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
0081     {
0082         using std::abs;
0083         // this overwrites x_err !
0084         algebra.for_each3( x_err , x_old , dxdt_old ,
0085                 typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * abs(get_unit_value( dt )) ) );
0086
0087         // value_type res = algebra.reduce( x_err ,
0088         //        typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0 ) );
0089         return algebra.norm_inf( x_err );
0090     }
0091
0092 private:
0093
0094     value_type m_eps_abs;
0095     value_type m_eps_rel;
0096     value_type m_a_x;
0097     value_type m_a_dxdt;
0098
0099 };
0100
0101
0102 template< typename Value, typename Time >
0103 class default_step_adjuster
0104 {
0105 public:
0106     typedef Time time_type;
0107     typedef Value value_type;
0108
0109     default_step_adjuster(const time_type max_dt=static_cast<time_type>(0))
0110             : m_max_dt(max_dt)
0111     {}
0112
0113
0114     time_type decrease_step(time_type dt, const value_type error, const int error_order) const
0115     {
0116         // returns the decreased time step
0117         BOOST_USING_STD_MIN();
0118         BOOST_USING_STD_MAX();
0119         using std::pow;
0120
0121         dt *= max
0122         BOOST_PREVENT_MACRO_SUBSTITUTION(
0123                 static_cast<value_type>( static_cast<value_type>(9) / static_cast<value_type>(10) *
0124                                          pow(error, static_cast<value_type>(-1) / (error_order - 1))),
0125                 static_cast<value_type>( static_cast<value_type>(1) / static_cast<value_type> (5)));
0126         if(m_max_dt != static_cast<time_type >(0))
0127             // limit to maximal stepsize even when decreasing
0128             dt = detail::min_abs(dt, m_max_dt);
0129         return dt;
0130     }
0131
0132     time_type increase_step(time_type dt, value_type error, const int stepper_order) const
0133     {
0134         // returns the increased time step
0135         BOOST_USING_STD_MIN();
0136         BOOST_USING_STD_MAX();
0137         using std::pow;
0138
0139         // adjust the size if dt is smaller than max_dt (providede max_dt is not zero)
0140         if(error < 0.5)
0141         {
0142             // error should be > 0
0143             error = max BOOST_PREVENT_MACRO_SUBSTITUTION (
0144                     static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(stepper_order) ) ) ,
0145                     error);
0146             // time_type dt_old = dt;   unused variable warning
0147             //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
0148             dt *= static_cast<value_type>(9)/static_cast<value_type>(10) *
0149                   pow(error, static_cast<value_type>(-1) / stepper_order);
0150             if(m_max_dt != static_cast<time_type >(0))
0151                 // limit to maximal stepsize
0152                 dt = detail::min_abs(dt, m_max_dt);
0153         }
0154         return dt;
0155     }
0156
0157     bool check_step_size_limit(const time_type dt)
0158     {
0159         if(m_max_dt != static_cast<time_type >(0))
0160             return detail::less_eq_with_sign(dt, m_max_dt, dt);
0161         return true;
0162     }
0163
0164     time_type get_max_dt() { return m_max_dt; }
0165
0166 protected:
0167     time_type m_max_dt;
0168 };
0169
0170
0171
0172 /*
0173  * error stepper category dispatcher
0174  */
0175 template<
0176 class ErrorStepper ,
0177 class ErrorChecker = default_error_checker< typename ErrorStepper::value_type ,
0178     typename ErrorStepper::algebra_type ,
0179     typename ErrorStepper::operations_type > ,
0180 class StepAdjuster = default_step_adjuster< typename ErrorStepper::value_type ,
0181     typename ErrorStepper::time_type > ,
0182 class Resizer = typename ErrorStepper::resizer_type ,
0183 class ErrorStepperCategory = typename ErrorStepper::stepper_category
0184 >
0185 class controlled_runge_kutta ;
0186
0187
0188
0189 /*
0190  * explicit stepper version
0191  *
0192  * this class introduces the following try_step overloads
0193     * try_step( sys , x , t , dt )
0194     * try_step( sys , x , dxdt , t , dt )
0195     * try_step( sys , in , t , out , dt )
0196     * try_step( sys , in , dxdt , t , out , dt )
0197  */
0198 /**
0199  * \brief Implements step size control for Runge-Kutta steppers with error
0200  * estimation.
0201  *
0202  * This class implements the step size control for standard Runge-Kutta
0203  * steppers with error estimation.
0204  *
0205  * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept.
0206  * \tparam ErrorChecker The error checker
0207  * \tparam Resizer The resizer policy type.
0208  */
0209 template<
0210 class ErrorStepper,
0211 class ErrorChecker,
0212 class StepAdjuster,
0213 class Resizer
0214 >
0215 class controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster, Resizer ,
0216         explicit_error_stepper_tag >
0217 {
0218
0219 public:
0220
0221     typedef ErrorStepper stepper_type;
0222     typedef typename stepper_type::state_type state_type;
0223     typedef typename stepper_type::value_type value_type;
0224     typedef typename stepper_type::deriv_type deriv_type;
0225     typedef typename stepper_type::time_type time_type;
0226     typedef typename stepper_type::algebra_type algebra_type;
0227     typedef typename stepper_type::operations_type operations_type;
0228     typedef Resizer resizer_type;
0229     typedef ErrorChecker error_checker_type;
0230     typedef StepAdjuster step_adjuster_type;
0231     typedef explicit_controlled_stepper_tag stepper_category;
0232
0233 #ifndef DOXYGEN_SKIP
0234     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
0235     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
0236
0237     typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster ,
0238             Resizer , explicit_error_stepper_tag > controlled_stepper_type;
0239 #endif //DOXYGEN_SKIP
0240
0241
0242     /**
0243      * \brief Constructs the controlled Runge-Kutta stepper.
0244      * \param error_checker An instance of the error checker.
0245      * \param stepper An instance of the underlying stepper.
0246      */
0247     controlled_runge_kutta(
0248             const error_checker_type &error_checker = error_checker_type( ) ,
0249             const step_adjuster_type &step_adjuster = step_adjuster_type() ,
0250             const stepper_type &stepper = stepper_type( )
0251     )
0252         : m_stepper(stepper), m_error_checker(error_checker) , m_step_adjuster(step_adjuster)
0253     { }
0254
0255
0256
0257     /*
0258      * Version 1 : try_step( sys , x , t , dt )
0259      *
0260      * The overloads are needed to solve the forwarding problem
0261      */
0262     /**
0263      * \brief Tries to perform one step.
0264      *
0265      * This method tries to do one step with step size dt. If the error estimate
0266      * is to large, the step is rejected and the method returns fail and the
0267      * step size dt is reduced. If the error estimate is acceptably small, the
0268      * step is performed, success is returned and dt might be increased to make
0269      * the steps as large as possible. This method also updates t if a step is
0270      * performed.
0271      *
0272      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0273      *               Simple System concept.
0274      * \param x The state of the ODE which should be solved. Overwritten if
0275      * the step is successful.
0276      * \param t The value of the time. Updated if the step is successful.
0277      * \param dt The step size. Updated.
0278      * \return success if the step was accepted, fail otherwise.
0279      */
0280     template< class System , class StateInOut >
0281     controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
0282     {
0283         return try_step_v1( system , x , t, dt );
0284     }
0285
0286     /**
0287      * \brief Tries to perform one step. Solves the forwarding problem and
0288      * allows for using boost range as state_type.
0289      *
0290      * This method tries to do one step with step size dt. If the error estimate
0291      * is to large, the step is rejected and the method returns fail and the
0292      * step size dt is reduced. If the error estimate is acceptably small, the
0293      * step is performed, success is returned and dt might be increased to make
0294      * the steps as large as possible. This method also updates t if a step is
0295      * performed.
0296      *
0297      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0298      *               Simple System concept.
0299      * \param x The state of the ODE which should be solved. Overwritten if
0300      * the step is successful. Can be a boost range.
0301      * \param t The value of the time. Updated if the step is successful.
0302      * \param dt The step size. Updated.
0303      * \return success if the step was accepted, fail otherwise.
0304      */
0305     template< class System , class StateInOut >
0306     controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
0307     {
0308         return try_step_v1( system , x , t, dt );
0309     }
0310
0311
0312
0313     /*
0314      * Version 2 : try_step( sys , x , dxdt , t , dt )
0315      *
0316      * this version does not solve the forwarding problem, boost.range can not be used
0317      */
0318     /**
0319      * \brief Tries to perform one step.
0320      *
0321      * This method tries to do one step with step size dt. If the error estimate
0322      * is to large, the step is rejected and the method returns fail and the
0323      * step size dt is reduced. If the error estimate is acceptably small, the
0324      * step is performed, success is returned and dt might be increased to make
0325      * the steps as large as possible. This method also updates t if a step is
0326      * performed.
0327      *
0328      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0329      *               Simple System concept.
0330      * \param x The state of the ODE which should be solved. Overwritten if
0331      * the step is successful.
0332      * \param dxdt The derivative of state.
0333      * \param t The value of the time. Updated if the step is successful.
0334      * \param dt The step size. Updated.
0335      * \return success if the step was accepted, fail otherwise.
0336      */
0337     template< class System , class StateInOut , class DerivIn >
0338     controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt )
0339     {
0340         m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
0341         controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt );
0342         if( res == success )
0343         {
0344             boost::numeric::odeint::copy( m_xnew.m_v , x );
0345         }
0346         return res;
0347     }
0348
0349     /*
0350      * Version 3 : try_step( sys , in , t , out , dt )
0351      *
0352      * this version does not solve the forwarding problem, boost.range can not be used
0353      *
0354      * the disable is needed to avoid ambiguous overloads if state_type = time_type
0355      */
0356     /**
0357      * \brief Tries to perform one step.
0358      *
0359      * \note This method is disabled if state_type=time_type to avoid ambiguity.
0360      *
0361      * This method tries to do one step with step size dt. If the error estimate
0362      * is to large, the step is rejected and the method returns fail and the
0363      * step size dt is reduced. If the error estimate is acceptably small, the
0364      * step is performed, success is returned and dt might be increased to make
0365      * the steps as large as possible. This method also updates t if a step is
0366      * performed.
0367      *
0368      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0369      *               Simple System concept.
0370      * \param in The state of the ODE which should be solved.
0371      * \param t The value of the time. Updated if the step is successful.
0372      * \param out Used to store the result of the step.
0373      * \param dt The step size. Updated.
0374      * \return success if the step was accepted, fail otherwise.
0375      */
0376     template< class System , class StateIn , class StateOut >
0377     typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type
0378     try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
0379     {
0380         typename odeint::unwrap_reference< System >::type &sys = system;
0381         m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
0382         sys( in , m_dxdt.m_v , t );
0383         return try_step( system , in , m_dxdt.m_v , t , out , dt );
0384     }
0385
0386
0387     /*
0388      * Version 4 : try_step( sys , in , dxdt , t , out , dt )
0389      *
0390      * this version does not solve the forwarding problem, boost.range can not be used
0391      */
0392     /**
0393      * \brief Tries to perform one step.
0394      *
0395      * This method tries to do one step with step size dt. If the error estimate
0396      * is to large, the step is rejected and the method returns fail and the
0397      * step size dt is reduced. If the error estimate is acceptably small, the
0398      * step is performed, success is returned and dt might be increased to make
0399      * the steps as large as possible. This method also updates t if a step is
0400      * performed.
0401      *
0402      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0403      *               Simple System concept.
0404      * \param in The state of the ODE which should be solved.
0405      * \param dxdt The derivative of state.
0406      * \param t The value of the time. Updated if the step is successful.
0407      * \param out Used to store the result of the step.
0408      * \param dt The step size. Updated.
0409      * \return success if the step was accepted, fail otherwise.
0410      */
0411     template< class System , class StateIn , class DerivIn , class StateOut >
0412     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
0413     {
0414         unwrapped_step_adjuster &step_adjuster = m_step_adjuster;
0415         if( !step_adjuster.check_step_size_limit(dt) )
0416         {
0417             // given dt was above step size limit - adjust and return fail;
0418             dt = step_adjuster.get_max_dt();
0419             return fail;
0420         }
0421
0422         m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
0423
0424         // do one step with error calculation
0425         m_stepper.do_step( system , in , dxdt , t , out , dt , m_xerr.m_v );
0426
0427         value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt , m_xerr.m_v , dt );
0428
0429         if( max_rel_err > 1.0 )
0430         {
0431             // error too big, decrease step size and reject this step
0432             dt = step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
0433             return fail;
0434         } else
0435         {
0436             // otherwise, increase step size and accept
0437             t += dt;
0438             dt = step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
0439             return success;
0440         }
0441     }
0442
0443     /**
0444      * \brief Adjust the size of all temporaries in the stepper manually.
0445      * \param x A state from which the size of the temporaries to be resized is deduced.
0446      */
0447     template< class StateType >
0448     void adjust_size( const StateType &x )
0449     {
0450         resize_m_xerr_impl( x );
0451         resize_m_dxdt_impl( x );
0452         resize_m_xnew_impl( x );
0453         m_stepper.adjust_size( x );
0454     }
0455
0456     /**
0457      * \brief Returns the instance of the underlying stepper.
0458      * \returns The instance of the underlying stepper.
0459      */
0460     stepper_type& stepper( void )
0461     {
0462         return m_stepper;
0463     }
0464
0465     /**
0466      * \brief Returns the instance of the underlying stepper.
0467      * \returns The instance of the underlying stepper.
0468      */
0469     const stepper_type& stepper( void ) const
0470     {
0471         return m_stepper;
0472     }
0473
0474 private:
0475
0476
0477     template< class System , class StateInOut >
0478     controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
0479     {
0480         typename odeint::unwrap_reference< System >::type &sys = system;
0481         m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
0482         sys( x , m_dxdt.m_v ,t );
0483         return try_step( system , x , m_dxdt.m_v , t , dt );
0484     }
0485
0486     template< class StateIn >
0487     bool resize_m_xerr_impl( const StateIn &x )
0488     {
0489         return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() );
0490     }
0491
0492     template< class StateIn >
0493     bool resize_m_dxdt_impl( const StateIn &x )
0494     {
0495         return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
0496     }
0497
0498     template< class StateIn >
0499     bool resize_m_xnew_impl( const StateIn &x )
0500     {
0501         return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
0502     }
0503
0504
0505
0506     stepper_type m_stepper;
0507     error_checker_type m_error_checker;
0508     step_adjuster_type m_step_adjuster;
0509     typedef typename unwrap_reference< step_adjuster_type >::type unwrapped_step_adjuster;
0510
0511     resizer_type m_dxdt_resizer;
0512     resizer_type m_xerr_resizer;
0513     resizer_type m_xnew_resizer;
0514
0515     wrapped_deriv_type m_dxdt;
0516     wrapped_state_type m_xerr;
0517     wrapped_state_type m_xnew;
0518 };
0519
0520
0521
0522
0523
0524
0525
0526
0527
0528
0529 /*
0530  * explicit stepper fsal version
0531  *
0532  * the class introduces the following try_step overloads
0533     * try_step( sys , x , t , dt )
0534     * try_step( sys , in , t , out , dt )
0535     * try_step( sys , x , dxdt , t , dt )
0536     * try_step( sys , in , dxdt_in , t , out , dxdt_out , dt )
0537  */
0538 /**
0539  * \brief Implements step size control for Runge-Kutta FSAL steppers with
0540  * error estimation.
0541  *
0542  * This class implements the step size control for FSAL Runge-Kutta
0543  * steppers with error estimation.
0544  *
0545  * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept.
0546  * \tparam ErrorChecker The error checker
0547  * \tparam Resizer The resizer policy type.
0548  */
0549 template<
0550 class ErrorStepper ,
0551 class ErrorChecker ,
0552 class StepAdjuster ,
0553 class Resizer
0554 >
0555 class controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster , Resizer , explicit_error_stepper_fsal_tag >
0556 {
0557
0558 public:
0559
0560     typedef ErrorStepper stepper_type;
0561     typedef typename stepper_type::state_type state_type;
0562     typedef typename stepper_type::value_type value_type;
0563     typedef typename stepper_type::deriv_type deriv_type;
0564     typedef typename stepper_type::time_type time_type;
0565     typedef typename stepper_type::algebra_type algebra_type;
0566     typedef typename stepper_type::operations_type operations_type;
0567     typedef Resizer resizer_type;
0568     typedef ErrorChecker error_checker_type;
0569     typedef StepAdjuster step_adjuster_type;
0570     typedef explicit_controlled_stepper_fsal_tag stepper_category;
0571
0572 #ifndef DOXYGEN_SKIP
0573     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
0574     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
0575
0576     typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
0577 #endif // DOXYGEN_SKIP
0578
0579     /**
0580      * \brief Constructs the controlled Runge-Kutta stepper.
0581      * \param error_checker An instance of the error checker.
0582      * \param stepper An instance of the underlying stepper.
0583      */
0584     controlled_runge_kutta(
0585             const error_checker_type &error_checker = error_checker_type() ,
0586             const step_adjuster_type &step_adjuster = step_adjuster_type() ,
0587             const stepper_type &stepper = stepper_type()
0588     )
0589     : m_stepper( stepper ) , m_error_checker( error_checker ) , m_step_adjuster(step_adjuster) ,
0590       m_first_call( true )
0591     { }
0592
0593     /*
0594      * Version 1 : try_step( sys , x , t , dt )
0595      *
0596      * The two overloads are needed in order to solve the forwarding problem
0597      */
0598     /**
0599      * \brief Tries to perform one step.
0600      *
0601      * This method tries to do one step with step size dt. If the error estimate
0602      * is to large, the step is rejected and the method returns fail and the
0603      * step size dt is reduced. If the error estimate is acceptably small, the
0604      * step is performed, success is returned and dt might be increased to make
0605      * the steps as large as possible. This method also updates t if a step is
0606      * performed.
0607      *
0608      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0609      *               Simple System concept.
0610      * \param x The state of the ODE which should be solved. Overwritten if
0611      * the step is successful.
0612      * \param t The value of the time. Updated if the step is successful.
0613      * \param dt The step size. Updated.
0614      * \return success if the step was accepted, fail otherwise.
0615      */
0616     template< class System , class StateInOut >
0617     controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
0618     {
0619         return try_step_v1( system , x , t , dt );
0620     }
0621
0622
0623     /**
0624      * \brief Tries to perform one step. Solves the forwarding problem and
0625      * allows for using boost range as state_type.
0626      *
0627      * This method tries to do one step with step size dt. If the error estimate
0628      * is to large, the step is rejected and the method returns fail and the
0629      * step size dt is reduced. If the error estimate is acceptably small, the
0630      * step is performed, success is returned and dt might be increased to make
0631      * the steps as large as possible. This method also updates t if a step is
0632      * performed.
0633      *
0634      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0635      *               Simple System concept.
0636      * \param x The state of the ODE which should be solved. Overwritten if
0637      * the step is successful. Can be a boost range.
0638      * \param t The value of the time. Updated if the step is successful.
0639      * \param dt The step size. Updated.
0640      * \return success if the step was accepted, fail otherwise.
0641      */
0642     template< class System , class StateInOut >
0643     controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
0644     {
0645         return try_step_v1( system , x , t , dt );
0646     }
0647
0648
0649
0650     /*
0651      * Version 2 : try_step( sys , in , t , out , dt );
0652      *
0653      * This version does not solve the forwarding problem, boost::range can not be used.
0654      *
0655      * The disabler is needed to solve ambiguous overloads
0656      */
0657     /**
0658      * \brief Tries to perform one step.
0659      *
0660      * \note This method is disabled if state_type=time_type to avoid ambiguity.
0661      *
0662      * This method tries to do one step with step size dt. If the error estimate
0663      * is to large, the step is rejected and the method returns fail and the
0664      * step size dt is reduced. If the error estimate is acceptably small, the
0665      * step is performed, success is returned and dt might be increased to make
0666      * the steps as large as possible. This method also updates t if a step is
0667      * performed.
0668      *
0669      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0670      *               Simple System concept.
0671      * \param in The state of the ODE which should be solved.
0672      * \param t The value of the time. Updated if the step is successful.
0673      * \param out Used to store the result of the step.
0674      * \param dt The step size. Updated.
0675      * \return success if the step was accepted, fail otherwise.
0676      */
0677     template< class System , class StateIn , class StateOut >
0678     typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type
0679     try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
0680     {
0681         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 )
0682         {
0683             initialize( system , in , t );
0684         }
0685         return try_step( system , in , m_dxdt.m_v , t , out , dt );
0686     }
0687
0688
0689     /*
0690      * Version 3 : try_step( sys , x , dxdt , t , dt )
0691      *
0692      * This version does not solve the forwarding problem, boost::range can not be used.
0693      */
0694     /**
0695      * \brief Tries to perform one step.
0696      *
0697      * This method tries to do one step with step size dt. If the error estimate
0698      * is to large, the step is rejected and the method returns fail and the
0699      * step size dt is reduced. If the error estimate is acceptably small, the
0700      * step is performed, success is returned and dt might be increased to make
0701      * the steps as large as possible. This method also updates t if a step is
0702      * performed.
0703      *
0704      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0705      *               Simple System concept.
0706      * \param x The state of the ODE which should be solved. Overwritten if
0707      * the step is successful.
0708      * \param dxdt The derivative of state.
0709      * \param t The value of the time. Updated if the step is successful.
0710      * \param dt The step size. Updated.
0711      * \return success if the step was accepted, fail otherwise.
0712      */
0713     template< class System , class StateInOut , class DerivInOut >
0714     controlled_step_result try_step( System system , StateInOut &x , DerivInOut &dxdt , time_type &t , time_type &dt )
0715     {
0716         m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
0717         m_dxdt_new_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_new_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
0718         controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , m_dxdtnew.m_v , dt );
0719         if( res == success )
0720         {
0721             boost::numeric::odeint::copy( m_xnew.m_v , x );
0722             boost::numeric::odeint::copy( m_dxdtnew.m_v , dxdt );
0723         }
0724         return res;
0725     }
0726
0727
0728     /*
0729      * Version 4 : try_step( sys , in , dxdt_in , t , out , dxdt_out , dt )
0730      *
0731      * This version does not solve the forwarding problem, boost::range can not be used.
0732      */
0733     /**
0734      * \brief Tries to perform one step.
0735      *
0736      * This method tries to do one step with step size dt. If the error estimate
0737      * is to large, the step is rejected and the method returns fail and the
0738      * step size dt is reduced. If the error estimate is acceptably small, the
0739      * step is performed, success is returned and dt might be increased to make
0740      * the steps as large as possible. This method also updates t if a step is
0741      * performed.
0742      *
0743      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0744      *               Simple System concept.
0745      * \param in The state of the ODE which should be solved.
0746      * \param dxdt The derivative of state.
0747      * \param t The value of the time. Updated if the step is successful.
0748      * \param out Used to store the result of the step.
0749      * \param dt The step size. Updated.
0750      * \return success if the step was accepted, fail otherwise.
0751      */
0752     template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
0753     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t ,
0754             StateOut &out , DerivOut &dxdt_out , time_type &dt )
0755     {
0756         unwrapped_step_adjuster &step_adjuster = m_step_adjuster;
0757         if( !step_adjuster.check_step_size_limit(dt) )
0758         {
0759             // given dt was above step size limit - adjust and return fail;
0760             dt = step_adjuster.get_max_dt();
0761             return fail;
0762         }
0763
0764         m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
0765
0766         //fsal: m_stepper.get_dxdt( dxdt );
0767         //fsal: m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err );
0768         m_stepper.do_step( system , in , dxdt_in , t , out , dxdt_out , dt , m_xerr.m_v );
0769
0770         // this potentially overwrites m_x_err! (standard_error_checker does, at least)
0771         value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt_in , m_xerr.m_v , dt );
0772
0773         if( max_rel_err > 1.0 )
0774         {
0775             // error too big, decrease step size and reject this step
0776             dt = step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
0777             return fail;
0778         }
0779         // otherwise, increase step size and accept
0780         t += dt;
0781         dt = step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
0782         return success;
0783     }
0784
0785
0786     /**
0787      * \brief Resets the internal state of the underlying FSAL stepper.
0788      */
0789     void reset( void )
0790     {
0791         m_first_call = true;
0792     }
0793
0794     /**
0795      * \brief Initializes the internal state storing an internal copy of the derivative.
0796      *
0797      * \param deriv The initial derivative of the ODE.
0798      */
0799     template< class DerivIn >
0800     void initialize( const DerivIn &deriv )
0801     {
0802         boost::numeric::odeint::copy( deriv , m_dxdt.m_v );
0803         m_first_call = false;
0804     }
0805
0806     /**
0807      * \brief Initializes the internal state storing an internal copy of the derivative.
0808      *
0809      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0810      *               Simple System concept.
0811      * \param x The initial state of the ODE which should be solved.
0812      * \param t The initial time.
0813      */
0814     template< class System , class StateIn >
0815     void initialize( System system , const StateIn &x , time_type t )
0816     {
0817         typename odeint::unwrap_reference< System >::type &sys = system;
0818         sys( x , m_dxdt.m_v , t );
0819         m_first_call = false;
0820     }
0821
0822     /**
0823      * \brief Returns true if the stepper has been initialized, false otherwise.
0824      *
0825      * \return true, if the stepper has been initialized, false otherwise.
0826      */
0827     bool is_initialized( void ) const
0828     {
0829         return ! m_first_call;
0830     }
0831
0832
0833     /**
0834      * \brief Adjust the size of all temporaries in the stepper manually.
0835      * \param x A state from which the size of the temporaries to be resized is deduced.
0836      */
0837     template< class StateType >
0838     void adjust_size( const StateType &x )
0839     {
0840         resize_m_xerr_impl( x );
0841         resize_m_dxdt_impl( x );
0842         resize_m_dxdt_new_impl( x );
0843         resize_m_xnew_impl( x );
0844     }
0845
0846
0847     /**
0848      * \brief Returns the instance of the underlying stepper.
0849      * \returns The instance of the underlying stepper.
0850      */
0851     stepper_type& stepper( void )
0852     {
0853         return m_stepper;
0854     }
0855
0856     /**
0857      * \brief Returns the instance of the underlying stepper.
0858      * \returns The instance of the underlying stepper.
0859      */
0860     const stepper_type& stepper( void ) const
0861     {
0862         return m_stepper;
0863     }
0864
0865
0866
0867 private:
0868
0869
0870     template< class StateIn >
0871     bool resize_m_xerr_impl( const StateIn &x )
0872     {
0873         return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() );
0874     }
0875
0876     template< class StateIn >
0877     bool resize_m_dxdt_impl( const StateIn &x )
0878     {
0879         return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
0880     }
0881
0882     template< class StateIn >
0883     bool resize_m_dxdt_new_impl( const StateIn &x )
0884     {
0885         return adjust_size_by_resizeability( m_dxdtnew , x , typename is_resizeable<deriv_type>::type() );
0886     }
0887
0888     template< class StateIn >
0889     bool resize_m_xnew_impl( const StateIn &x )
0890     {
0891         return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
0892     }
0893
0894
0895     template< class System , class StateInOut >
0896     controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
0897     {
0898         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 )
0899         {
0900             initialize( system , x , t );
0901         }
0902         return try_step( system , x , m_dxdt.m_v , t , dt );
0903     }
0904
0905
0906     stepper_type m_stepper;
0907     error_checker_type m_error_checker;
0908     step_adjuster_type m_step_adjuster;
0909     typedef typename unwrap_reference< step_adjuster_type >::type unwrapped_step_adjuster;
0910
0911     resizer_type m_dxdt_resizer;
0912     resizer_type m_xerr_resizer;
0913     resizer_type m_xnew_resizer;
0914     resizer_type m_dxdt_new_resizer;
0915
0916     wrapped_deriv_type m_dxdt;
0917     wrapped_state_type m_xerr;
0918     wrapped_state_type m_xnew;
0919     wrapped_deriv_type m_dxdtnew;
0920     bool m_first_call;
0921 };
0922
0923
0924 /********** DOXYGEN **********/
0925
0926 /**** DEFAULT ERROR CHECKER ****/
0927
0928 /**
0929  * \class default_error_checker
0930  * \brief The default error checker to be used with Runge-Kutta error steppers
0931  *
0932  * This class provides the default mechanism to compare the error estimates
0933  * reported by Runge-Kutta error steppers with user defined error bounds.
0934  * It is used by the controlled_runge_kutta steppers.
0935  *
0936  * \tparam Value The value type.
0937  * \tparam Time The time type.
0938  * \tparam Algebra The algebra type.
0939  * \tparam Operations The operations type.
0940  */
0941
0942     /**
0943      * \fn default_error_checker( value_type eps_abs , value_type eps_rel , value_type a_x , value_type a_dxdt ,
0944      * time_type max_dt)
0945      * \brief Constructs the error checker.
0946      *
0947      * The error is calculated as follows: ????
0948      *
0949      * \param eps_abs Absolute tolerance level.
0950      * \param eps_rel Relative tolerance level.
0951      * \param a_x Factor for the weight of the state.
0952      * \param a_dxdt Factor for the weight of the derivative.
0953      * \param max_dt Maximum allowed step size.
0954      */
0955
0956     /**
0957      * \fn error( const State &x_old , const Deriv &dxdt_old , Err &x_err , time_type dt ) const
0958      * \brief Calculates the error level.
0959      *
0960      * If the returned error level is greater than 1, the estimated error was
0961      * larger than the permitted error bounds and the step should be repeated
0962      * with a smaller step size.
0963      *
0964      * \param x_old State at the beginning of the step.
0965      * \param dxdt_old Derivative at the beginning of the step.
0966      * \param x_err Error estimate.
0967      * \param dt Time step.
0968      * \return error
0969      */
0970
0971     /**
0972      * \fn error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , time_type dt ) const
0973      * \brief Calculates the error level using a given algebra.
0974      *
0975      * If the returned error level is greater than 1, the estimated error was
0976      * larger than the permitted error bounds and the step should be repeated
0977      * with a smaller step size.
0978      *
0979      * \param algebra The algebra used for calculation of the error.
0980      * \param x_old State at the beginning of the step.
0981      * \param dxdt_old Derivative at the beginning of the step.
0982      * \param x_err Error estimate.
0983      * \param dt Time step.
0984      * \return error
0985      */
0986
0987     /**
0988      * \fn time_type decrease_step(const time_type dt, const value_type error, const int error_order)
0989      * \brief Returns a decreased step size based on the given error and order
0990      *
0991      * Calculates a new smaller step size based on the given error and its order.
0992      *
0993      * \param dt The old step size.
0994      * \param error The computed error estimate.
0995      * \param error_order The error order of the stepper.
0996      * \return dt_new The new, reduced step size.
0997      */
0998
0999     /**
1000      * \fn time_type increase_step(const time_type dt, const value_type error, const int error_order)
1001      * \brief Returns an increased step size based on the given error and order.
1002      *
1003      * Calculates a new bigger step size based on the given error and its order. If max_dt != 0, the
1004      * new step size is limited to max_dt.
1005      *
1006      * \param dt The old step size.
1007      * \param error The computed error estimate.
1008      * \param error_order The order of the stepper.
1009      * \return dt_new The new, increased step size.
1010      */
1011
1012
1013 } // odeint
1014 } // numeric
1015 } // boost
1016
1017
1018 #endif // BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED