ceph/src/boost/libs/numeric/odeint/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp

   1 /* [auto_generated]
   2  boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
   3
   4  [begin_description]
   5  The default controlled stepper which can be used with all explicit Runge-Kutta error steppers.
   6  [end_description]
   7
   8  Copyright 2010-2013 Karsten Ahnert
   9  Copyright 2010-2015 Mario Mulansky
  10  Copyright 2012 Christoph Koke
  11
  12  Distributed under the Boost Software License, Version 1.0.
  13  (See accompanying file LICENSE_1_0.txt or
  14  copy at http://www.boost.org/LICENSE_1_0.txt)
  15  */
  16
  17
  18 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED
  19 #define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED
  20
  21
  22
  23 #include <cmath>
  24
  25 #include <boost/config.hpp>
  26 #include <boost/utility/enable_if.hpp>
  27 #include <boost/type_traits/is_same.hpp>
  28
  29 #include <boost/numeric/odeint/util/bind.hpp>
  30 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
  31 #include <boost/numeric/odeint/util/copy.hpp>
  32
  33 #include <boost/numeric/odeint/util/state_wrapper.hpp>
  34 #include <boost/numeric/odeint/util/is_resizeable.hpp>
  35 #include <boost/numeric/odeint/util/resizer.hpp>
  36 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
  37
  38 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
  39 #include <boost/numeric/odeint/algebra/default_operations.hpp>
  40 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
  41
  42 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
  43 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
  44
  45 namespace boost {
  46 namespace numeric {
  47 namespace odeint {
  48
  49
  50 template
  51 <
  52 class Value ,
  53 class Algebra ,
  54 class Operations
  55 >
  56 class default_error_checker
  57 {
  58 public:
  59
  60     typedef Value value_type;
  61     typedef Algebra algebra_type;
  62     typedef Operations operations_type;
  63
  64     default_error_checker(
  65             value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
  66             value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
  67             value_type a_x = static_cast< value_type >( 1 ) ,
  68             value_type a_dxdt = static_cast< value_type >( 1 ))
  69         : m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
  70     { }
  71
  72
  73     template< class State , class Deriv , class Err, class Time >
  74     value_type error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
  75     {
  76         return error( algebra_type() , x_old , dxdt_old , x_err , dt );
  77     }
  78
  79     template< class State , class Deriv , class Err, class Time >
  80     value_type error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
  81     {
  82         using std::abs;
  83         // this overwrites x_err !
  84         algebra.for_each3( x_err , x_old , dxdt_old ,
  85                 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 )) ) );
  86
  87         // value_type res = algebra.reduce( x_err ,
  88         //        typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0 ) );
  89         return algebra.norm_inf( x_err );
  90     }
  91
  92 private:
  93
  94     value_type m_eps_abs;
  95     value_type m_eps_rel;
  96     value_type m_a_x;
  97     value_type m_a_dxdt;
  98
  99 };
 100
 101
 102 template< typename Value, typename Time >
 103 class default_step_adjuster
 104 {
 105 public:
 106     typedef Time time_type;
 107     typedef Value value_type;
 108
 109     default_step_adjuster(const time_type max_dt=static_cast<time_type>(0))
 110             : m_max_dt(max_dt)
 111     {}
 112
 113
 114     time_type decrease_step(time_type dt, const value_type error, const int error_order) const
 115     {
 116         // returns the decreased time step
 117         BOOST_USING_STD_MIN();
 118         BOOST_USING_STD_MAX();
 119         using std::pow;
 120
 121         dt *= max
 122         BOOST_PREVENT_MACRO_SUBSTITUTION(
 123                 static_cast<value_type>( static_cast<value_type>(9) / static_cast<value_type>(10) *
 124                                          pow(error, static_cast<value_type>(-1) / (error_order - 1))),
 125                 static_cast<value_type>( static_cast<value_type>(1) / static_cast<value_type> (5)));
 126         if(m_max_dt != static_cast<time_type >(0))
 127             // limit to maximal stepsize even when decreasing
 128             dt = detail::min_abs(dt, m_max_dt);
 129         return dt;
 130     }
 131
 132     time_type increase_step(time_type dt, value_type error, const int stepper_order) const
 133     {
 134         // returns the increased time step
 135         BOOST_USING_STD_MIN();
 136         BOOST_USING_STD_MAX();
 137         using std::pow;
 138
 139         // adjust the size if dt is smaller than max_dt (providede max_dt is not zero)
 140         if(error < 0.5)
 141         {
 142             // error should be > 0
 143             error = max BOOST_PREVENT_MACRO_SUBSTITUTION (
 144                     static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(stepper_order) ) ) ,
 145                     error);
 146             time_type dt_old = dt;
 147             //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
 148             dt *= static_cast<value_type>(9)/static_cast<value_type>(10) *
 149                   pow(error, static_cast<value_type>(-1) / stepper_order);
 150             if(m_max_dt != static_cast<time_type >(0))
 151                 // limit to maximal stepsize
 152                 dt = detail::min_abs(dt, m_max_dt);
 153         }
 154         return dt;
 155     }
 156
 157     bool check_step_size_limit(const time_type dt)
 158     {
 159         if(m_max_dt != static_cast<time_type >(0))
 160             return detail::less_eq_with_sign(dt, m_max_dt, dt);
 161         return true;
 162     }
 163
 164     time_type get_max_dt() { return m_max_dt; }
 165
 166 private:
 167     time_type m_max_dt;
 168 };
 169
 170
 171
 172 /*
 173  * error stepper category dispatcher
 174  */
 175 template<
 176 class ErrorStepper ,
 177 class ErrorChecker = default_error_checker< typename ErrorStepper::value_type ,
 178     typename ErrorStepper::algebra_type ,
 179     typename ErrorStepper::operations_type > ,
 180 class StepAdjuster = default_step_adjuster< typename ErrorStepper::value_type ,
 181     typename ErrorStepper::time_type > ,
 182 class Resizer = typename ErrorStepper::resizer_type ,
 183 class ErrorStepperCategory = typename ErrorStepper::stepper_category
 184 >
 185 class controlled_runge_kutta ;
 186
 187
 188
 189 /*
 190  * explicit stepper version
 191  *
 192  * this class introduces the following try_step overloads
 193     * try_step( sys , x , t , dt )
 194     * try_step( sys , x , dxdt , t , dt )
 195     * try_step( sys , in , t , out , dt )
 196     * try_step( sys , in , dxdt , t , out , dt )
 197  */
 198 /**
 199  * \brief Implements step size control for Runge-Kutta steppers with error
 200  * estimation.
 201  *
 202  * This class implements the step size control for standard Runge-Kutta
 203  * steppers with error estimation.
 204  *
 205  * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept.
 206  * \tparam ErrorChecker The error checker
 207  * \tparam Resizer The resizer policy type.
 208  */
 209 template<
 210 class ErrorStepper,
 211 class ErrorChecker,
 212 class StepAdjuster,
 213 class Resizer
 214 >
 215 class controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster, Resizer ,
 216         explicit_error_stepper_tag >
 217 {
 218
 219 public:
 220
 221     typedef ErrorStepper stepper_type;
 222     typedef typename stepper_type::state_type state_type;
 223     typedef typename stepper_type::value_type value_type;
 224     typedef typename stepper_type::deriv_type deriv_type;
 225     typedef typename stepper_type::time_type time_type;
 226     typedef typename stepper_type::algebra_type algebra_type;
 227     typedef typename stepper_type::operations_type operations_type;
 228     typedef Resizer resizer_type;
 229     typedef ErrorChecker error_checker_type;
 230     typedef StepAdjuster step_adjuster_type;
 231     typedef explicit_controlled_stepper_tag stepper_category;
 232
 233 #ifndef DOXYGEN_SKIP
 234     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
 235     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
 236
 237     typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster ,
 238             Resizer , explicit_error_stepper_tag > controlled_stepper_type;
 239 #endif //DOXYGEN_SKIP
 240
 241
 242     /**
 243      * \brief Constructs the controlled Runge-Kutta stepper.
 244      * \param error_checker An instance of the error checker.
 245      * \param stepper An instance of the underlying stepper.
 246      */
 247     controlled_runge_kutta(
 248             const error_checker_type &error_checker = error_checker_type( ) ,
 249             const step_adjuster_type &step_adjuster = step_adjuster_type() ,
 250             const stepper_type &stepper = stepper_type( )
 251     )
 252         : m_stepper(stepper), m_error_checker(error_checker) , m_step_adjuster(step_adjuster)
 253     { }
 254
 255
 256
 257     /*
 258      * Version 1 : try_step( sys , x , t , dt )
 259      *
 260      * The overloads are needed to solve the forwarding problem
 261      */
 262     /**
 263      * \brief Tries to perform one step.
 264      *
 265      * This method tries to do one step with step size dt. If the error estimate
 266      * is to large, the step is rejected and the method returns fail and the
 267      * step size dt is reduced. If the error estimate is acceptably small, the
 268      * step is performed, success is returned and dt might be increased to make
 269      * the steps as large as possible. This method also updates t if a step is
 270      * performed.
 271      *
 272      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
 273      *               Simple System concept.
 274      * \param x The state of the ODE which should be solved. Overwritten if
 275      * the step is successful.
 276      * \param t The value of the time. Updated if the step is successful.
 277      * \param dt The step size. Updated.
 278      * \return success if the step was accepted, fail otherwise.
 279      */
 280     template< class System , class StateInOut >
 281     controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
 282     {
 283         return try_step_v1( system , x , t, dt );
 284     }
 285
 286     /**
 287      * \brief Tries to perform one step. Solves the forwarding problem and
 288      * allows for using boost range as state_type.
 289      *
 290      * This method tries to do one step with step size dt. If the error estimate
 291      * is to large, the step is rejected and the method returns fail and the
 292      * step size dt is reduced. If the error estimate is acceptably small, the
 293      * step is performed, success is returned and dt might be increased to make
 294      * the steps as large as possible. This method also updates t if a step is
 295      * performed.
 296      *
 297      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
 298      *               Simple System concept.
 299      * \param x The state of the ODE which should be solved. Overwritten if
 300      * the step is successful. Can be a boost range.
 301      * \param t The value of the time. Updated if the step is successful.
 302      * \param dt The step size. Updated.
 303      * \return success if the step was accepted, fail otherwise.
 304      */
 305     template< class System , class StateInOut >
 306     controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
 307     {
 308         return try_step_v1( system , x , t, dt );
 309     }
 310
 311
 312
 313     /*
 314      * Version 2 : try_step( sys , x , dxdt , t , dt )
 315      *
 316      * this version does not solve the forwarding problem, boost.range can not be used
 317      */
 318     /**
 319      * \brief Tries to perform one step.
 320      *
 321      * This method tries to do one step with step size dt. If the error estimate
 322      * is to large, the step is rejected and the method returns fail and the
 323      * step size dt is reduced. If the error estimate is acceptably small, the
 324      * step is performed, success is returned and dt might be increased to make
 325      * the steps as large as possible. This method also updates t if a step is
 326      * performed.
 327      *
 328      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
 329      *               Simple System concept.
 330      * \param x The state of the ODE which should be solved. Overwritten if
 331      * the step is successful.
 332      * \param dxdt The derivative of state.
 333      * \param t The value of the time. Updated if the step is successful.
 334      * \param dt The step size. Updated.
 335      * \return success if the step was accepted, fail otherwise.
 336      */
 337     template< class System , class StateInOut , class DerivIn >
 338     controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt )
 339     {
 340         m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
 341         controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt );
 342         if( res == success )
 343         {
 344             boost::numeric::odeint::copy( m_xnew.m_v , x );
 345         }
 346         return res;
 347     }
 348
 349     /*
 350      * Version 3 : try_step( sys , in , t , out , dt )
 351      *
 352      * this version does not solve the forwarding problem, boost.range can not be used
 353      *
 354      * the disable is needed to avoid ambiguous overloads if state_type = time_type
 355      */
 356     /**
 357      * \brief Tries to perform one step.
 358      *
 359      * \note This method is disabled if state_type=time_type to avoid ambiguity.
 360      *
 361      * This method tries to do one step with step size dt. If the error estimate
 362      * is to large, the step is rejected and the method returns fail and the
 363      * step size dt is reduced. If the error estimate is acceptably small, the
 364      * step is performed, success is returned and dt might be increased to make
 365      * the steps as large as possible. This method also updates t if a step is
 366      * performed.
 367      *
 368      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
 369      *               Simple System concept.
 370      * \param in The state of the ODE which should be solved.
 371      * \param t The value of the time. Updated if the step is successful.
 372      * \param out Used to store the result of the step.
 373      * \param dt The step size. Updated.
 374      * \return success if the step was accepted, fail otherwise.
 375      */
 376     template< class System , class StateIn , class StateOut >
 377     typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type
 378     try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
 379     {
 380         typename odeint::unwrap_reference< System >::type &sys = system;
 381         m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
 382         sys( in , m_dxdt.m_v , t );
 383         return try_step( system , in , m_dxdt.m_v , t , out , dt );
 384     }
 385
 386
 387     /*
 388      * Version 4 : try_step( sys , in , dxdt , t , out , dt )
 389      *
 390      * this version does not solve the forwarding problem, boost.range can not be used
 391      */
 392     /**
 393      * \brief Tries to perform one step.
 394      *
 395      * This method tries to do one step with step size dt. If the error estimate
 396      * is to large, the step is rejected and the method returns fail and the
 397      * step size dt is reduced. If the error estimate is acceptably small, the
 398      * step is performed, success is returned and dt might be increased to make
 399      * the steps as large as possible. This method also updates t if a step is
 400      * performed.
 401      *
 402      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
 403      *               Simple System concept.
 404      * \param in The state of the ODE which should be solved.
 405      * \param dxdt The derivative of state.
 406      * \param t The value of the time. Updated if the step is successful.
 407      * \param out Used to store the result of the step.
 408      * \param dt The step size. Updated.
 409      * \return success if the step was accepted, fail otherwise.
 410      */
 411     template< class System , class StateIn , class DerivIn , class StateOut >
 412     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
 413     {
 414         if( !m_step_adjuster.check_step_size_limit(dt) )
 415         {
 416             // given dt was above step size limit - adjust and return fail;
 417             dt = m_step_adjuster.get_max_dt();
 418             return fail;
 419         }
 420
 421         m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
 422
 423         // do one step with error calculation
 424         m_stepper.do_step( system , in , dxdt , t , out , dt , m_xerr.m_v );
 425
 426         value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt , m_xerr.m_v , dt );
 427
 428         if( max_rel_err > 1.0 )
 429         {
 430             // error too big, decrease step size and reject this step
 431             dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
 432             return fail;
 433         } else
 434         {
 435             // otherwise, increase step size and accept
 436             t += dt;
 437             dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
 438             return success;
 439         }
 440     }
 441
 442     /**
 443      * \brief Adjust the size of all temporaries in the stepper manually.
 444      * \param x A state from which the size of the temporaries to be resized is deduced.
 445      */
 446     template< class StateType >
 447     void adjust_size( const StateType &x )
 448     {
 449         resize_m_xerr_impl( x );
 450         resize_m_dxdt_impl( x );
 451         resize_m_xnew_impl( x );
 452         m_stepper.adjust_size( x );
 453     }
 454
 455     /**
 456      * \brief Returns the instance of the underlying stepper.
 457      * \returns The instance of the underlying stepper.
 458      */
 459     stepper_type& stepper( void )
 460     {
 461         return m_stepper;
 462     }
 463
 464     /**
 465      * \brief Returns the instance of the underlying stepper.
 466      * \returns The instance of the underlying stepper.
 467      */
 468     const stepper_type& stepper( void ) const
 469     {
 470         return m_stepper;
 471     }
 472
 473 private:
 474
 475
 476     template< class System , class StateInOut >
 477     controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
 478     {
 479         typename odeint::unwrap_reference< System >::type &sys = system;
 480         m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
 481         sys( x , m_dxdt.m_v ,t );
 482         return try_step( system , x , m_dxdt.m_v , t , dt );
 483     }
 484
 485     template< class StateIn >
 486     bool resize_m_xerr_impl( const StateIn &x )
 487     {
 488         return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() );
 489     }
 490
 491     template< class StateIn >
 492     bool resize_m_dxdt_impl( const StateIn &x )
 493     {
 494         return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
 495     }
 496
 497     template< class StateIn >
 498     bool resize_m_xnew_impl( const StateIn &x )
 499     {
 500         return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
 501     }
 502
 503
 504
 505     stepper_type m_stepper;
 506     error_checker_type m_error_checker;
 507     step_adjuster_type m_step_adjuster;
 508
 509     resizer_type m_dxdt_resizer;
 510     resizer_type m_xerr_resizer;
 511     resizer_type m_xnew_resizer;
 512
 513     wrapped_deriv_type m_dxdt;
 514     wrapped_state_type m_xerr;
 515     wrapped_state_type m_xnew;
 516 };
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527 /*
 528  * explicit stepper fsal version
 529  *
 530  * the class introduces the following try_step overloads
 531     * try_step( sys , x , t , dt )
 532     * try_step( sys , in , t , out , dt )
 533     * try_step( sys , x , dxdt , t , dt )
 534     * try_step( sys , in , dxdt_in , t , out , dxdt_out , dt )
 535  */
 536 /**
 537  * \brief Implements step size control for Runge-Kutta FSAL steppers with
 538  * error estimation.
 539  *
 540  * This class implements the step size control for FSAL Runge-Kutta
 541  * steppers with error estimation.
 542  *
 543  * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept.
 544  * \tparam ErrorChecker The error checker
 545  * \tparam Resizer The resizer policy type.
 546  */
 547 template<
 548 class ErrorStepper ,
 549 class ErrorChecker ,
 550 class StepAdjuster ,
 551 class Resizer
 552 >
 553 class controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster , Resizer , explicit_error_stepper_fsal_tag >
 554 {
 555
 556 public:
 557
 558     typedef ErrorStepper stepper_type;
 559     typedef typename stepper_type::state_type state_type;
 560     typedef typename stepper_type::value_type value_type;
 561     typedef typename stepper_type::deriv_type deriv_type;
 562     typedef typename stepper_type::time_type time_type;
 563     typedef typename stepper_type::algebra_type algebra_type;
 564     typedef typename stepper_type::operations_type operations_type;
 565     typedef Resizer resizer_type;
 566     typedef ErrorChecker error_checker_type;
 567     typedef StepAdjuster step_adjuster_type;
 568     typedef explicit_controlled_stepper_fsal_tag stepper_category;
 569
 570 #ifndef DOXYGEN_SKIP
 571     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
 572     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
 573
 574     typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
 575 #endif // DOXYGEN_SKIP
 576
 577     /**
 578      * \brief Constructs the controlled Runge-Kutta stepper.
 579      * \param error_checker An instance of the error checker.
 580      * \param stepper An instance of the underlying stepper.
 581      */
 582     controlled_runge_kutta(
 583             const error_checker_type &error_checker = error_checker_type() ,
 584             const step_adjuster_type &step_adjuster = step_adjuster_type() ,
 585             const stepper_type &stepper = stepper_type()
 586     )
 587     : m_stepper( stepper ) , m_error_checker( error_checker ) , m_step_adjuster(step_adjuster) ,
 588       m_first_call( true )
 589     { }
 590
 591     /*
 592      * Version 1 : try_step( sys , x , t , dt )
 593      *
 594      * The two overloads are needed in order to solve the forwarding problem
 595      */
 596     /**
 597      * \brief Tries to perform one step.
 598      *
 599      * This method tries to do one step with step size dt. If the error estimate
 600      * is to large, the step is rejected and the method returns fail and the
 601      * step size dt is reduced. If the error estimate is acceptably small, the
 602      * step is performed, success is returned and dt might be increased to make
 603      * the steps as large as possible. This method also updates t if a step is
 604      * performed.
 605      *
 606      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
 607      *               Simple System concept.
 608      * \param x The state of the ODE which should be solved. Overwritten if
 609      * the step is successful.
 610      * \param t The value of the time. Updated if the step is successful.
 611      * \param dt The step size. Updated.
 612      * \return success if the step was accepted, fail otherwise.
 613      */
 614     template< class System , class StateInOut >
 615     controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
 616     {
 617         return try_step_v1( system , x , t , dt );
 618     }
 619
 620
 621     /**
 622      * \brief Tries to perform one step. Solves the forwarding problem and
 623      * allows for using boost range as state_type.
 624      *
 625      * This method tries to do one step with step size dt. If the error estimate
 626      * is to large, the step is rejected and the method returns fail and the
 627      * step size dt is reduced. If the error estimate is acceptably small, the
 628      * step is performed, success is returned and dt might be increased to make
 629      * the steps as large as possible. This method also updates t if a step is
 630      * performed.
 631      *
 632      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
 633      *               Simple System concept.
 634      * \param x The state of the ODE which should be solved. Overwritten if
 635      * the step is successful. Can be a boost range.
 636      * \param t The value of the time. Updated if the step is successful.
 637      * \param dt The step size. Updated.
 638      * \return success if the step was accepted, fail otherwise.
 639      */
 640     template< class System , class StateInOut >
 641     controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
 642     {
 643         return try_step_v1( system , x , t , dt );
 644     }
 645
 646
 647
 648     /*
 649      * Version 2 : try_step( sys , in , t , out , dt );
 650      *
 651      * This version does not solve the forwarding problem, boost::range can not be used.
 652      *
 653      * The disabler is needed to solve ambiguous overloads
 654      */
 655     /**
 656      * \brief Tries to perform one step.
 657      *
 658      * \note This method is disabled if state_type=time_type to avoid ambiguity.
 659      *
 660      * This method tries to do one step with step size dt. If the error estimate
 661      * is to large, the step is rejected and the method returns fail and the
 662      * step size dt is reduced. If the error estimate is acceptably small, the
 663      * step is performed, success is returned and dt might be increased to make
 664      * the steps as large as possible. This method also updates t if a step is
 665      * performed.
 666      *
 667      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
 668      *               Simple System concept.
 669      * \param in The state of the ODE which should be solved.
 670      * \param t The value of the time. Updated if the step is successful.
 671      * \param out Used to store the result of the step.
 672      * \param dt The step size. Updated.
 673      * \return success if the step was accepted, fail otherwise.
 674      */
 675     template< class System , class StateIn , class StateOut >
 676     typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type
 677     try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
 678     {
 679         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 )
 680         {
 681             initialize( system , in , t );
 682         }
 683         return try_step( system , in , m_dxdt.m_v , t , out , dt );
 684     }
 685
 686
 687     /*
 688      * Version 3 : try_step( sys , x , dxdt , t , dt )
 689      *
 690      * This version does not solve the forwarding problem, boost::range can not be used.
 691      */
 692     /**
 693      * \brief Tries to perform one step.
 694      *
 695      * This method tries to do one step with step size dt. If the error estimate
 696      * is to large, the step is rejected and the method returns fail and the
 697      * step size dt is reduced. If the error estimate is acceptably small, the
 698      * step is performed, success is returned and dt might be increased to make
 699      * the steps as large as possible. This method also updates t if a step is
 700      * performed.
 701      *
 702      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
 703      *               Simple System concept.
 704      * \param x The state of the ODE which should be solved. Overwritten if
 705      * the step is successful.
 706      * \param dxdt The derivative of state.
 707      * \param t The value of the time. Updated if the step is successful.
 708      * \param dt The step size. Updated.
 709      * \return success if the step was accepted, fail otherwise.
 710      */
 711     template< class System , class StateInOut , class DerivInOut >
 712     controlled_step_result try_step( System system , StateInOut &x , DerivInOut &dxdt , time_type &t , time_type &dt )
 713     {
 714         m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
 715         m_dxdt_new_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_new_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
 716         controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , m_dxdtnew.m_v , dt );
 717         if( res == success )
 718         {
 719             boost::numeric::odeint::copy( m_xnew.m_v , x );
 720             boost::numeric::odeint::copy( m_dxdtnew.m_v , dxdt );
 721         }
 722         return res;
 723     }
 724
 725
 726     /*
 727      * Version 4 : try_step( sys , in , dxdt_in , t , out , dxdt_out , dt )
 728      *
 729      * This version does not solve the forwarding problem, boost::range can not be used.
 730      */
 731     /**
 732      * \brief Tries to perform one step.
 733      *
 734      * This method tries to do one step with step size dt. If the error estimate
 735      * is to large, the step is rejected and the method returns fail and the
 736      * step size dt is reduced. If the error estimate is acceptably small, the
 737      * step is performed, success is returned and dt might be increased to make
 738      * the steps as large as possible. This method also updates t if a step is
 739      * performed.
 740      *
 741      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
 742      *               Simple System concept.
 743      * \param in The state of the ODE which should be solved.
 744      * \param dxdt The derivative of state.
 745      * \param t The value of the time. Updated if the step is successful.
 746      * \param out Used to store the result of the step.
 747      * \param dt The step size. Updated.
 748      * \return success if the step was accepted, fail otherwise.
 749      */
 750     template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
 751     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t ,
 752             StateOut &out , DerivOut &dxdt_out , time_type &dt )
 753     {
 754         if( !m_step_adjuster.check_step_size_limit(dt) )
 755         {
 756             // given dt was above step size limit - adjust and return fail;
 757             dt = m_step_adjuster.get_max_dt();
 758             return fail;
 759         }
 760
 761         m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
 762
 763         //fsal: m_stepper.get_dxdt( dxdt );
 764         //fsal: m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err );
 765         m_stepper.do_step( system , in , dxdt_in , t , out , dxdt_out , dt , m_xerr.m_v );
 766
 767         // this potentially overwrites m_x_err! (standard_error_checker does, at least)
 768         value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt_in , m_xerr.m_v , dt );
 769
 770         if( max_rel_err > 1.0 )
 771         {
 772             // error too big, decrease step size and reject this step
 773             dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
 774             return fail;
 775         }
 776         // otherwise, increase step size and accept
 777         t += dt;
 778         dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
 779         return success;
 780     }
 781
 782
 783     /**
 784      * \brief Resets the internal state of the underlying FSAL stepper.
 785      */
 786     void reset( void )
 787     {
 788         m_first_call = true;
 789     }
 790
 791     /**
 792      * \brief Initializes the internal state storing an internal copy of the derivative.
 793      *
 794      * \param deriv The initial derivative of the ODE.
 795      */
 796     template< class DerivIn >
 797     void initialize( const DerivIn &deriv )
 798     {
 799         boost::numeric::odeint::copy( deriv , m_dxdt.m_v );
 800         m_first_call = false;
 801     }
 802
 803     /**
 804      * \brief Initializes the internal state storing an internal copy of the derivative.
 805      *
 806      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
 807      *               Simple System concept.
 808      * \param x The initial state of the ODE which should be solved.
 809      * \param t The initial time.
 810      */
 811     template< class System , class StateIn >
 812     void initialize( System system , const StateIn &x , time_type t )
 813     {
 814         typename odeint::unwrap_reference< System >::type &sys = system;
 815         sys( x , m_dxdt.m_v , t );
 816         m_first_call = false;
 817     }
 818
 819     /**
 820      * \brief Returns true if the stepper has been initialized, false otherwise.
 821      *
 822      * \return true, if the stepper has been initialized, false otherwise.
 823      */
 824     bool is_initialized( void ) const
 825     {
 826         return ! m_first_call;
 827     }
 828
 829
 830     /**
 831      * \brief Adjust the size of all temporaries in the stepper manually.
 832      * \param x A state from which the size of the temporaries to be resized is deduced.
 833      */
 834     template< class StateType >
 835     void adjust_size( const StateType &x )
 836     {
 837         resize_m_xerr_impl( x );
 838         resize_m_dxdt_impl( x );
 839         resize_m_dxdt_new_impl( x );
 840         resize_m_xnew_impl( x );
 841     }
 842
 843
 844     /**
 845      * \brief Returns the instance of the underlying stepper.
 846      * \returns The instance of the underlying stepper.
 847      */
 848     stepper_type& stepper( void )
 849     {
 850         return m_stepper;
 851     }
 852
 853     /**
 854      * \brief Returns the instance of the underlying stepper.
 855      * \returns The instance of the underlying stepper.
 856      */
 857     const stepper_type& stepper( void ) const
 858     {
 859         return m_stepper;
 860     }
 861
 862
 863
 864 private:
 865
 866
 867     template< class StateIn >
 868     bool resize_m_xerr_impl( const StateIn &x )
 869     {
 870         return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() );
 871     }
 872
 873     template< class StateIn >
 874     bool resize_m_dxdt_impl( const StateIn &x )
 875     {
 876         return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
 877     }
 878
 879     template< class StateIn >
 880     bool resize_m_dxdt_new_impl( const StateIn &x )
 881     {
 882         return adjust_size_by_resizeability( m_dxdtnew , x , typename is_resizeable<deriv_type>::type() );
 883     }
 884
 885     template< class StateIn >
 886     bool resize_m_xnew_impl( const StateIn &x )
 887     {
 888         return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
 889     }
 890
 891
 892     template< class System , class StateInOut >
 893     controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
 894     {
 895         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 )
 896         {
 897             initialize( system , x , t );
 898         }
 899         return try_step( system , x , m_dxdt.m_v , t , dt );
 900     }
 901
 902
 903     stepper_type m_stepper;
 904     error_checker_type m_error_checker;
 905     step_adjuster_type m_step_adjuster;
 906
 907     resizer_type m_dxdt_resizer;
 908     resizer_type m_xerr_resizer;
 909     resizer_type m_xnew_resizer;
 910     resizer_type m_dxdt_new_resizer;
 911
 912     wrapped_deriv_type m_dxdt;
 913     wrapped_state_type m_xerr;
 914     wrapped_state_type m_xnew;
 915     wrapped_deriv_type m_dxdtnew;
 916     bool m_first_call;
 917 };
 918
 919
 920 /********** DOXYGEN **********/
 921
 922 /**** DEFAULT ERROR CHECKER ****/
 923
 924 /**
 925  * \class default_error_checker
 926  * \brief The default error checker to be used with Runge-Kutta error steppers
 927  *
 928  * This class provides the default mechanism to compare the error estimates
 929  * reported by Runge-Kutta error steppers with user defined error bounds.
 930  * It is used by the controlled_runge_kutta steppers.
 931  *
 932  * \tparam Value The value type.
 933  * \tparam Time The time type.
 934  * \tparam Algebra The algebra type.
 935  * \tparam Operations The operations type.
 936  */
 937
 938     /**
 939      * \fn default_error_checker( value_type eps_abs , value_type eps_rel , value_type a_x , value_type a_dxdt ,
 940      * time_type max_dt)
 941      * \brief Constructs the error checker.
 942      *
 943      * The error is calculated as follows: ????
 944      *
 945      * \param eps_abs Absolute tolerance level.
 946      * \param eps_rel Relative tolerance level.
 947      * \param a_x Factor for the weight of the state.
 948      * \param a_dxdt Factor for the weight of the derivative.
 949      * \param max_dt Maximum allowed step size.
 950      */
 951
 952     /**
 953      * \fn error( const State &x_old , const Deriv &dxdt_old , Err &x_err , time_type dt ) const
 954      * \brief Calculates the error level.
 955      *
 956      * If the returned error level is greater than 1, the estimated error was
 957      * larger than the permitted error bounds and the step should be repeated
 958      * with a smaller step size.
 959      *
 960      * \param x_old State at the beginning of the step.
 961      * \param dxdt_old Derivative at the beginning of the step.
 962      * \param x_err Error estimate.
 963      * \param dt Time step.
 964      * \return error
 965      */
 966
 967     /**
 968      * \fn error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , time_type dt ) const
 969      * \brief Calculates the error level using a given algebra.
 970      *
 971      * If the returned error level is greater than 1, the estimated error was
 972      * larger than the permitted error bounds and the step should be repeated
 973      * with a smaller step size.
 974      *
 975      * \param algebra The algebra used for calculation of the error.
 976      * \param x_old State at the beginning of the step.
 977      * \param dxdt_old Derivative at the beginning of the step.
 978      * \param x_err Error estimate.
 979      * \param dt Time step.
 980      * \return error
 981      */
 982
 983     /**
 984      * \fn time_type decrease_step(const time_type dt, const value_type error, const int error_order)
 985      * \brief Returns a decreased step size based on the given error and order
 986      *
 987      * Calculates a new smaller step size based on the given error and its order.
 988      *
 989      * \param dt The old step size.
 990      * \param error The computed error estimate.
 991      * \param error_order The error order of the stepper.
 992      * \return dt_new The new, reduced step size.
 993      */
 994
 995     /**
 996      * \fn time_type increase_step(const time_type dt, const value_type error, const int error_order)
 997      * \brief Returns an increased step size based on the given error and order.
 998      *
 999      * Calculates a new bigger step size based on the given error and its order. If max_dt != 0, the
1000      * new step size is limited to max_dt.
1001      *
1002      * \param dt The old step size.
1003      * \param error The computed error estimate.
1004      * \param error_order The order of the stepper.
1005      * \return dt_new The new, increased step size.
1006      */
1007
1008
1009 } // odeint
1010 } // numeric
1011 } // boost
1012
1013
1014 #endif // BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED