ceph/src/boost/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp

   1 /*
   2  [auto_generated]
   3  boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp
   4
   5  [begin_description]
   6  Implementation of the Dense-output stepper for all steppers. Note, that this class does
   7  not computes the result but serves as an interface.
   8  [end_description]
   9
  10  Copyright 2011-2013 Karsten Ahnert
  11  Copyright 2011-2015 Mario Mulansky
  12  Copyright 2012 Christoph Koke
  13
  14  Distributed under the Boost Software License, Version 1.0.
  15  (See accompanying file LICENSE_1_0.txt or
  16  copy at http://www.boost.org/LICENSE_1_0.txt)
  17  */
  18
  19
  20 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED
  21 #define BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED
  22
  23
  24 #include <utility>
  25 #include <stdexcept>
  26
  27 #include <boost/throw_exception.hpp>
  28
  29 #include <boost/numeric/odeint/util/bind.hpp>
  30
  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
  37 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
  38 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
  39
  40 #include <boost/numeric/odeint/integrate/max_step_checker.hpp>
  41
  42 namespace boost {
  43 namespace numeric {
  44 namespace odeint {
  45
  46 template< class Stepper , class StepperCategory = typename Stepper::stepper_category >
  47 class dense_output_runge_kutta;
  48
  49
  50 /**
  51  * \brief The class representing dense-output Runge-Kutta steppers.
  52  * \note In this stepper, the initialize method has to be called before using
  53  * the do_step method.
  54  *
  55  * The dense-output functionality allows to interpolate the solution between
  56  * subsequent integration points using intermediate results obtained during the
  57  * computation. This version works based on a normal stepper without step-size
  58  * control.
  59  *
  60  *
  61  * \tparam Stepper The stepper type of the underlying algorithm.
  62  */
  63 template< class Stepper >
  64 class dense_output_runge_kutta< Stepper , stepper_tag >
  65 {
  66
  67 public:
  68
  69     /*
  70      * We do not need all typedefs.
  71      */
  72     typedef Stepper stepper_type;
  73     typedef typename stepper_type::state_type state_type;
  74     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
  75     typedef typename stepper_type::value_type value_type;
  76     typedef typename stepper_type::deriv_type deriv_type;
  77     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
  78     typedef typename stepper_type::time_type time_type;
  79     typedef typename stepper_type::algebra_type algebra_type;
  80     typedef typename stepper_type::operations_type operations_type;
  81     typedef typename stepper_type::resizer_type resizer_type;
  82     typedef dense_output_stepper_tag stepper_category;
  83     typedef dense_output_runge_kutta< Stepper > dense_output_stepper_type;
  84
  85
  86     /**
  87      * \brief Constructs the dense_output_runge_kutta class. An instance of the
  88      * underlying stepper can be provided.
  89      * \param stepper An instance of the underlying stepper.
  90      */
  91     dense_output_runge_kutta( const stepper_type &stepper = stepper_type() )
  92     : m_stepper( stepper ) , m_resizer() ,
  93       m_x1() , m_x2() , m_current_state_x1( true ) ,
  94       m_t() , m_t_old() , m_dt()
  95     { }
  96
  97
  98     /**
  99      * \brief Initializes the stepper. Has to be called before do_step can be
 100      * used to set the initial conditions and the step size.
 101      * \param x0 The initial state of the ODE which should be solved.
 102      * \param t0 The initial time, at which the step should be performed.
 103      * \param dt0 The step size.
 104      */
 105     template< class StateType >
 106     void initialize( const StateType &x0 , time_type t0 , time_type dt0 )
 107     {
 108         m_resizer.adjust_size( x0 , detail::bind( &dense_output_stepper_type::template resize_impl< StateType > , detail::ref( *this ) , detail::_1 ) );
 109         boost::numeric::odeint::copy( x0 , get_current_state() );
 110         m_t = t0;
 111         m_dt = dt0;
 112     }
 113
 114     /**
 115      * \brief Does one time step.
 116      * \note initialize has to be called before using this method to set the
 117      * initial conditions x,t and the stepsize.
 118      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
 119      *               Simple System concept.
 120      * \return Pair with start and end time of the integration step.
 121      */
 122     template< class System >
 123     std::pair< time_type , time_type > do_step( System system )
 124     {
 125         m_stepper.do_step( system , get_current_state() , m_t , get_old_state() , m_dt );
 126         m_t_old = m_t;
 127         m_t += m_dt;
 128         toggle_current_state();
 129         return std::make_pair( m_t_old , m_t );
 130     }
 131
 132     /*
 133      * The next two overloads are needed to solve the forwarding problem
 134      */
 135
 136     /**
 137      * \brief Calculates the solution at an intermediate point.
 138      * \param t The time at which the solution should be calculated, has to be
 139      * in the current time interval.
 140      * \param x The output variable where the result is written into.
 141      */
 142     template< class StateOut >
 143     void calc_state( time_type t , StateOut &x ) const
 144     {
 145         if( t == current_time() )
 146         {
 147             boost::numeric::odeint::copy( get_current_state() , x );
 148         }
 149         m_stepper.calc_state( x , t , get_old_state() , m_t_old , get_current_state() , m_t );
 150     }
 151
 152     /**
 153      * \brief Calculates the solution at an intermediate point. Solves the forwarding problem
 154      * \param t The time at which the solution should be calculated, has to be
 155      * in the current time interval.
 156      * \param x The output variable where the result is written into, can be a boost range.
 157      */
 158     template< class StateOut >
 159     void calc_state( time_type t , const StateOut &x ) const
 160     {
 161         m_stepper.calc_state( x , t , get_old_state() , m_t_old , get_current_state() , m_t );
 162     }
 163
 164     /**
 165      * \brief Adjust the size of all temporaries in the stepper manually.
 166      * \param x A state from which the size of the temporaries to be resized is deduced.
 167      */
 168     template< class StateType >
 169     void adjust_size( const StateType &x )
 170     {
 171         resize_impl( x );
 172         m_stepper.stepper().resize( x );
 173     }
 174
 175     /**
 176      * \brief Returns the current state of the solution.
 177      * \return The current state of the solution x(t).
 178      */
 179     const state_type& current_state( void ) const
 180     {
 181         return get_current_state();
 182     }
 183
 184     /**
 185      * \brief Returns the current time of the solution.
 186      * \return The current time of the solution t.
 187      */
 188     time_type current_time( void ) const
 189     {
 190         return m_t;
 191     }
 192
 193     /**
 194      * \brief Returns the last state of the solution.
 195      * \return The last state of the solution x(t-dt).
 196      */
 197     const state_type& previous_state( void ) const
 198     {
 199         return get_old_state();
 200     }
 201
 202     /**
 203      * \brief Returns the last time of the solution.
 204      * \return The last time of the solution t-dt.
 205      */
 206     time_type previous_time( void ) const
 207     {
 208         return m_t_old;
 209     }
 210
 211     /**
 212      * \brief Returns the current time step.
 213      * \return dt.
 214      */
 215     time_type current_time_step( void ) const
 216     {
 217         return m_dt;
 218     }
 219
 220
 221 private:
 222
 223     state_type& get_current_state( void )
 224     {
 225         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
 226     }
 227
 228     const state_type& get_current_state( void ) const
 229     {
 230         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
 231     }
 232
 233     state_type& get_old_state( void )
 234     {
 235         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
 236     }
 237
 238     const state_type& get_old_state( void ) const
 239     {
 240         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
 241     }
 242
 243     void toggle_current_state( void )
 244     {
 245         m_current_state_x1 = ! m_current_state_x1;
 246     }
 247
 248
 249     template< class StateIn >
 250     bool resize_impl( const StateIn &x )
 251     {
 252         bool resized = false;
 253         resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() );
 254         resized |= adjust_size_by_resizeability( m_x2 , x , typename is_resizeable<state_type>::type() );
 255         return resized;
 256     }
 257
 258
 259     stepper_type m_stepper;
 260     resizer_type m_resizer;
 261     wrapped_state_type m_x1 , m_x2;
 262     bool m_current_state_x1;    // if true, the current state is m_x1
 263     time_type m_t , m_t_old , m_dt;
 264
 265 };
 266
 267
 268
 269
 270
 271 /**
 272  * \brief The class representing dense-output Runge-Kutta steppers with FSAL property.
 273  *
 274  * The interface is the same as for dense_output_runge_kutta< Stepper , stepper_tag >.
 275  * This class provides dense output functionality based on methods with step size controlled
 276  *
 277  *
 278  * \tparam Stepper The stepper type of the underlying algorithm.
 279  */
 280 template< class Stepper >
 281 class dense_output_runge_kutta< Stepper , explicit_controlled_stepper_fsal_tag >
 282 {
 283 public:
 284
 285     /*
 286      * We do not need all typedefs.
 287      */
 288     typedef Stepper controlled_stepper_type;
 289
 290     typedef typename controlled_stepper_type::stepper_type stepper_type;
 291     typedef typename stepper_type::state_type state_type;
 292     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
 293     typedef typename stepper_type::value_type value_type;
 294     typedef typename stepper_type::deriv_type deriv_type;
 295     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
 296     typedef typename stepper_type::time_type time_type;
 297     typedef typename stepper_type::algebra_type algebra_type;
 298     typedef typename stepper_type::operations_type operations_type;
 299     typedef typename stepper_type::resizer_type resizer_type;
 300     typedef dense_output_stepper_tag stepper_category;
 301     typedef dense_output_runge_kutta< Stepper > dense_output_stepper_type;
 302
 303
 304     dense_output_runge_kutta( const controlled_stepper_type &stepper = controlled_stepper_type() )
 305     : m_stepper( stepper ) , m_resizer() ,
 306       m_current_state_x1( true ) ,
 307       m_x1() , m_x2() , m_dxdt1() , m_dxdt2() ,
 308       m_t() , m_t_old() , m_dt() ,
 309       m_is_deriv_initialized( false )
 310     { }
 311
 312
 313     template< class StateType >
 314     void initialize( const StateType &x0 , time_type t0 , time_type dt0 )
 315     {
 316         m_resizer.adjust_size( x0 , detail::bind( &dense_output_stepper_type::template resize< StateType > , detail::ref( *this ) , detail::_1 ) );
 317         boost::numeric::odeint::copy( x0 , get_current_state() );
 318         m_t = t0;
 319         m_dt = dt0;
 320         m_is_deriv_initialized = false;
 321     }
 322
 323     template< class System >
 324     std::pair< time_type , time_type > do_step( System system )
 325     {
 326         if( !m_is_deriv_initialized )
 327         {
 328             typename odeint::unwrap_reference< System >::type &sys = system;
 329             sys( get_current_state() , get_current_deriv() , m_t );
 330             m_is_deriv_initialized = true;
 331         }
 332
 333         failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
 334         controlled_step_result res = fail;
 335         m_t_old = m_t;
 336         do
 337         {
 338             res = m_stepper.try_step( system , get_current_state() , get_current_deriv() , m_t ,
 339                                       get_old_state() , get_old_deriv() , m_dt );
 340             fail_checker();  // check for overflow of failed steps
 341         }
 342         while( res == fail );
 343         toggle_current_state();
 344         return std::make_pair( m_t_old , m_t );
 345     }
 346
 347
 348     /*
 349      * The two overloads are needed in order to solve the forwarding problem.
 350      */
 351     template< class StateOut >
 352     void calc_state( time_type t , StateOut &x ) const
 353     {
 354         m_stepper.stepper().calc_state( t , x , get_old_state() , get_old_deriv() , m_t_old ,
 355                                         get_current_state() , get_current_deriv() , m_t );
 356     }
 357
 358     template< class StateOut >
 359     void calc_state( time_type t , const StateOut &x ) const
 360     {
 361         m_stepper.stepper().calc_state( t , x , get_old_state() , get_old_deriv() , m_t_old ,
 362                                         get_current_state() , get_current_deriv() , m_t );
 363     }
 364
 365
 366     template< class StateIn >
 367     bool resize( const StateIn &x )
 368     {
 369         bool resized = false;
 370         resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() );
 371         resized |= adjust_size_by_resizeability( m_x2 , x , typename is_resizeable<state_type>::type() );
 372         resized |= adjust_size_by_resizeability( m_dxdt1 , x , typename is_resizeable<deriv_type>::type() );
 373         resized |= adjust_size_by_resizeability( m_dxdt2 , x , typename is_resizeable<deriv_type>::type() );
 374         return resized;
 375     }
 376
 377
 378     template< class StateType >
 379     void adjust_size( const StateType &x )
 380     {
 381         resize( x );
 382         m_stepper.stepper().resize( x );
 383     }
 384
 385     const state_type& current_state( void ) const
 386     {
 387         return get_current_state();
 388     }
 389
 390     time_type current_time( void ) const
 391     {
 392         return m_t;
 393     }
 394
 395     const state_type& previous_state( void ) const
 396     {
 397         return get_old_state();
 398     }
 399
 400     time_type previous_time( void ) const
 401     {
 402         return m_t_old;
 403     }
 404
 405     time_type current_time_step( void ) const
 406     {
 407         return m_dt;
 408     }
 409
 410
 411 private:
 412
 413     state_type& get_current_state( void )
 414     {
 415         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
 416     }
 417
 418     const state_type& get_current_state( void ) const
 419     {
 420         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
 421     }
 422
 423     state_type& get_old_state( void )
 424     {
 425         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
 426     }
 427
 428     const state_type& get_old_state( void ) const
 429     {
 430         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
 431     }
 432
 433     deriv_type& get_current_deriv( void )
 434     {
 435         return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
 436     }
 437
 438     const deriv_type& get_current_deriv( void ) const
 439     {
 440         return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
 441     }
 442
 443     deriv_type& get_old_deriv( void )
 444     {
 445         return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
 446     }
 447
 448     const deriv_type& get_old_deriv( void ) const
 449     {
 450         return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
 451     }
 452
 453
 454     void toggle_current_state( void )
 455     {
 456         m_current_state_x1 = ! m_current_state_x1;
 457     }
 458
 459
 460     controlled_stepper_type m_stepper;
 461     resizer_type m_resizer;
 462     bool m_current_state_x1;
 463     wrapped_state_type m_x1 , m_x2;
 464     wrapped_deriv_type m_dxdt1 , m_dxdt2;
 465     time_type m_t , m_t_old , m_dt;
 466     bool m_is_deriv_initialized;
 467
 468 };
 469
 470 } // namespace odeint
 471 } // namespace numeric
 472 } // namespace boost
 473
 474
 475
 476 #endif // BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED