odeint/stepper/adams_bashforth.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/stepper/adams_bashforth.hpp
0004
0005  [begin_description]
0006  Implementaton of the Adam-Bashforth method a multistep method used for the predictor step in the
0007  Adams-Bashforth-Moulton method.
0008  [end_description]
0009
0010  Copyright 2011-2013 Karsten Ahnert
0011  Copyright 2011-2013 Mario Mulansky
0012  Copyright 2012 Christoph Koke
0013  Copyright 2013 Pascal Germroth
0014
0015  Distributed under the Boost Software License, Version 1.0.
0016  (See accompanying file LICENSE_1_0.txt or
0017  copy at http://www.boost.org/LICENSE_1_0.txt)
0018  */
0019
0020
0021 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_ADAMS_BASHFORTH_HPP_INCLUDED
0022 #define BOOST_NUMERIC_ODEINT_STEPPER_ADAMS_BASHFORTH_HPP_INCLUDED
0023
0024 #include <boost/static_assert.hpp>
0025
0026 #include <boost/numeric/odeint/util/bind.hpp>
0027 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0028
0029 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
0030 #include <boost/numeric/odeint/algebra/default_operations.hpp>
0031 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
0032 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
0033
0034 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0035 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0036 #include <boost/numeric/odeint/util/resizer.hpp>
0037
0038 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0039 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
0040 #include <boost/numeric/odeint/stepper/extrapolation_stepper.hpp>
0041
0042 #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
0043
0044 #include <boost/numeric/odeint/stepper/detail/adams_bashforth_coefficients.hpp>
0045 #include <boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp>
0046 #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
0047
0048 #include <boost/mpl/arithmetic.hpp>
0049 #include <boost/mpl/min_max.hpp>
0050 #include <boost/mpl/equal_to.hpp>
0051
0052 namespace mpl = boost::mpl;
0053
0054
0055 namespace boost {
0056 namespace numeric {
0057 namespace odeint {
0058
0059     using mpl::int_;
0060
0061     /* if N >= 4, returns the smallest even number > N, otherwise returns 4 */
0062     template < int N >
0063     struct order_helper
0064         : mpl::max< typename mpl::eval_if<
0065                         mpl::equal_to< mpl::modulus< int_< N >, int_< 2 > >,
0066                                        int_< 0 > >,
0067                         int_< N >, int_< N + 1 > >::type,
0068                     int_< 4 > >::type
0069     { };
0070
0071 template<
0072 size_t Steps ,
0073 class State ,
0074 class Value = double ,
0075 class Deriv = State ,
0076 class Time = Value ,
0077 class Algebra = typename algebra_dispatcher< State >::algebra_type ,
0078 class Operations = typename operations_dispatcher< State >::operations_type ,
0079 class Resizer = initially_resizer ,
0080 class InitializingStepper = extrapolation_stepper< order_helper<Steps>::value,
0081                                                    State, Value, Deriv, Time,
0082                                                    Algebra, Operations, Resizer >
0083 >
0084 class adams_bashforth : public algebra_stepper_base< Algebra , Operations >
0085 {
0086
0087 #ifndef DOXYGEN_SKIP
0088     BOOST_STATIC_ASSERT(( Steps > 0 ));
0089     BOOST_STATIC_ASSERT(( Steps < 9 ));
0090 #endif
0091
0092 public :
0093
0094     typedef State state_type;
0095     typedef state_wrapper< state_type > wrapped_state_type;
0096     typedef Value value_type;
0097     typedef Deriv deriv_type;
0098     typedef state_wrapper< deriv_type > wrapped_deriv_type;
0099     typedef Time time_type;
0100     typedef Resizer resizer_type;
0101     typedef stepper_tag stepper_category;
0102
0103     typedef InitializingStepper initializing_stepper_type;
0104
0105     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
0106     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
0107     typedef typename algebra_stepper_base_type::operations_type operations_type;
0108 #ifndef DOXYGEN_SKIP
0109     typedef adams_bashforth< Steps , State , Value , Deriv , Time , Algebra , Operations , Resizer , InitializingStepper > stepper_type;
0110 #endif
0111     static const size_t steps = Steps;
0112
0113
0114
0115     typedef unsigned short order_type;
0116     static const order_type order_value = steps;
0117
0118     typedef detail::rotating_buffer< wrapped_deriv_type , steps > step_storage_type;
0119
0120
0121
0122     order_type order( void ) const { return order_value; }
0123
0124     adams_bashforth( const algebra_type &algebra = algebra_type() )
0125     : algebra_stepper_base_type( algebra ) ,
0126       m_step_storage() , m_resizer() , m_coefficients() ,
0127       m_steps_initialized( 0 ) , m_initializing_stepper()
0128     { }
0129
0130
0131
0132     /*
0133      * Version 1 : do_step( system , x , t , dt );
0134      *
0135      * solves the forwarding problem
0136      */
0137     template< class System , class StateInOut >
0138     void do_step( System system , StateInOut &x , time_type t , time_type dt )
0139     {
0140         do_step( system , x , t , x , dt );
0141     }
0142
0143     /**
0144      * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateInOut.
0145      */
0146     template< class System , class StateInOut >
0147     void do_step( System system , const StateInOut &x , time_type t , time_type dt )
0148     {
0149         do_step( system , x , t , x , dt );
0150     }
0151
0152
0153
0154     /*
0155      * Version 2 : do_step( system , in , t , out , dt );
0156      *
0157      * solves the forwarding problem
0158      */
0159
0160     template< class System , class StateIn , class StateOut >
0161     void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
0162     {
0163         do_step_impl( system , in , t , out , dt );
0164     }
0165
0166     /**
0167      * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateOut.
0168      */
0169     template< class System , class StateIn , class StateOut >
0170     void do_step( System system , const StateIn &in , time_type t , const StateOut &out , time_type dt )
0171     {
0172         do_step_impl( system , in , t , out , dt );
0173     }
0174
0175
0176     template< class StateType >
0177     void adjust_size( const StateType &x )
0178     {
0179         resize_impl( x );
0180     }
0181
0182     const step_storage_type& step_storage( void ) const
0183     {
0184         return m_step_storage;
0185     }
0186
0187     step_storage_type& step_storage( void )
0188     {
0189         return m_step_storage;
0190     }
0191
0192     template< class ExplicitStepper , class System , class StateIn >
0193     void initialize( ExplicitStepper explicit_stepper , System system , StateIn &x , time_type &t , time_type dt )
0194     {
0195         typename odeint::unwrap_reference< ExplicitStepper >::type &stepper = explicit_stepper;
0196         typename odeint::unwrap_reference< System >::type &sys = system;
0197
0198         m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) );
0199
0200         for( size_t i=0 ; i+1<steps ; ++i )
0201         {
0202             if( i != 0 ) m_step_storage.rotate();
0203             sys( x , m_step_storage[0].m_v , t );
0204             stepper.do_step_dxdt_impl( system, x, m_step_storage[0].m_v, t,
0205                                        dt );
0206             t += dt;
0207         }
0208         m_steps_initialized = steps;
0209     }
0210
0211     template< class System , class StateIn >
0212     void initialize( System system , StateIn &x , time_type &t , time_type dt )
0213     {
0214         initialize( detail::ref( m_initializing_stepper ) , system , x , t , dt );
0215     }
0216
0217     void reset( void )
0218     {
0219         m_steps_initialized = 0;
0220     }
0221
0222     bool is_initialized( void ) const
0223     {
0224         return m_steps_initialized >= ( steps - 1 );
0225     }
0226
0227     const initializing_stepper_type& initializing_stepper( void ) const { return m_initializing_stepper; }
0228
0229     initializing_stepper_type& initializing_stepper( void ) { return m_initializing_stepper; }
0230
0231 private:
0232
0233     template< class System , class StateIn , class StateOut >
0234     void do_step_impl( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
0235     {
0236         typename odeint::unwrap_reference< System >::type &sys = system;
0237         if( m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) ) )
0238         {
0239             m_steps_initialized = 0;
0240         }
0241
0242         if( m_steps_initialized + 1 < steps )
0243         {
0244             if( m_steps_initialized != 0 ) m_step_storage.rotate();
0245             sys( in , m_step_storage[0].m_v , t );
0246             m_initializing_stepper.do_step_dxdt_impl(
0247                 system, in, m_step_storage[0].m_v, t, out, dt );
0248             ++m_steps_initialized;
0249         }
0250         else
0251         {
0252             m_step_storage.rotate();
0253             sys( in , m_step_storage[0].m_v , t );
0254             detail::adams_bashforth_call_algebra< steps , algebra_type , operations_type >()( this->m_algebra , in , out , m_step_storage , m_coefficients , dt );
0255         }
0256     }
0257
0258
0259     template< class StateIn >
0260     bool resize_impl( const StateIn &x )
0261     {
0262         bool resized( false );
0263         for( size_t i=0 ; i<steps ; ++i )
0264         {
0265             resized |= adjust_size_by_resizeability( m_step_storage[i] , x , typename is_resizeable<deriv_type>::type() );
0266         }
0267         return resized;
0268     }
0269
0270     step_storage_type m_step_storage;
0271     resizer_type m_resizer;
0272     detail::adams_bashforth_coefficients< value_type , steps > m_coefficients;
0273     size_t m_steps_initialized;
0274     initializing_stepper_type m_initializing_stepper;
0275
0276 };
0277
0278
0279 /***** DOXYGEN *****/
0280
0281 /**
0282  * \class adams_bashforth
0283  * \brief The Adams-Bashforth multistep algorithm.
0284  *
0285  * The Adams-Bashforth method is a multi-step algorithm with configurable step
0286  * number. The step number is specified as template parameter Steps and it
0287  * then uses the result from the previous Steps steps. See also
0288  * <a href="http://en.wikipedia.org/wiki/Linear_multistep_method">en.wikipedia.org/wiki/Linear_multistep_method</a>.
0289  * Currently, a maximum of Steps=8 is supported.
0290  * The method is explicit and fulfills the Stepper concept. Step size control
0291  * or continuous output are not provided.
0292  *
0293  * This class derives from algebra_base and inherits its interface via
0294  * CRTP (current recurring template pattern). For more details see
0295  * algebra_stepper_base.
0296  *
0297  * \tparam Steps The number of steps (maximal 8).
0298  * \tparam State The state type.
0299  * \tparam Value The value type.
0300  * \tparam Deriv The type representing the time derivative of the state.
0301  * \tparam Time The time representing the independent variable - the time.
0302  * \tparam Algebra The algebra type.
0303  * \tparam Operations The operations type.
0304  * \tparam Resizer The resizer policy type.
0305  * \tparam InitializingStepper The stepper for the first two steps.
0306  */
0307
0308     /**
0309      * \fn adams_bashforth::adams_bashforth( const algebra_type &algebra )
0310      * \brief Constructs the adams_bashforth class. This constructor can be used as a default
0311      * constructor if the algebra has a default constructor.
0312      * \param algebra A copy of algebra is made and stored.
0313      */
0314
0315     /**
0316      * \fn order_type adams_bashforth::order( void ) const
0317      * \brief Returns the order of the algorithm, which is equal to the number of steps.
0318      * \return order of the method.
0319      */
0320
0321     /**
0322      * \fn void adams_bashforth::do_step( System system , StateInOut &x , time_type t , time_type dt )
0323      * \brief This method performs one step. It transforms the result in-place.
0324      *
0325      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
0326      *               Simple System concept.
0327      * \param x The state of the ODE which should be solved. After calling do_step the result is updated in x.
0328      * \param t The value of the time, at which the step should be performed.
0329      * \param dt The step size.
0330      */
0331
0332     /**
0333      * \fn void adams_bashforth::do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
0334      * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
0335      *
0336      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
0337      *               Simple System concept.
0338      * \param in The state of the ODE which should be solved. in is not modified in this method
0339      * \param t The value of the time, at which the step should be performed.
0340      * \param out The result of the step is written in out.
0341      * \param dt The step size.
0342      */
0343
0344     /**
0345      * \fn void adams_bashforth::adjust_size( const StateType &x )
0346      * \brief Adjust the size of all temporaries in the stepper manually.
0347      * \param x A state from which the size of the temporaries to be resized is deduced.
0348      */
0349
0350
0351     /**
0352      * \fn const step_storage_type& adams_bashforth::step_storage( void ) const
0353      * \brief Returns the storage of intermediate results.
0354      * \return The storage of intermediate results.
0355      */
0356
0357     /**
0358      * \fn step_storage_type& adams_bashforth::step_storage( void )
0359      * \brief Returns the storage of intermediate results.
0360      * \return The storage of intermediate results.
0361      */
0362
0363     /**
0364      * \fn void adams_bashforth::initialize( ExplicitStepper explicit_stepper , System system , StateIn &x , time_type &t , time_type dt )
0365      * \brief Initialized the stepper. Does Steps-1 steps with the explicit_stepper to fill the buffer.
0366      * \param explicit_stepper the stepper used to fill the buffer of previous step results
0367      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
0368      *               Simple System concept.
0369      * \param x The state of the ODE which should be solved. After calling do_step the result is updated in x.
0370      * \param t The value of the time, at which the step should be performed.
0371      * \param dt The step size.
0372      */
0373
0374     /**
0375      * \fn void adams_bashforth::initialize( System system , StateIn &x , time_type &t , time_type dt )
0376      * \brief Initialized the stepper. Does Steps-1 steps with an internal instance of InitializingStepper to fill the buffer.
0377      * \note The state x and time t are updated to the values after Steps-1 initial steps.
0378      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
0379      *               Simple System concept.
0380      * \param x The initial state of the ODE which should be solved, updated in this method.
0381      * \param t The initial value of the time, updated in this method.
0382      * \param dt The step size.
0383      */
0384
0385     /**
0386      * \fn void adams_bashforth::reset( void )
0387      * \brief Resets the internal buffer of the stepper.
0388      */
0389
0390     /**
0391      * \fn bool adams_bashforth::is_initialized( void ) const
0392      * \brief Returns true if the stepper has been initialized.
0393      * \return bool true if stepper is initialized, false otherwise
0394      */
0395
0396     /**
0397      * \fn const initializing_stepper_type& adams_bashforth::initializing_stepper( void ) const
0398      * \brief Returns the internal initializing stepper instance.
0399      * \return initializing_stepper
0400      */
0401
0402     /**
0403      * \fn const initializing_stepper_type& adams_bashforth::initializing_stepper( void ) const
0404      * \brief Returns the internal initializing stepper instance.
0405      * \return initializing_stepper
0406      */
0407
0408     /**
0409      * \fn initializing_stepper_type& adams_bashforth::initializing_stepper( void )
0410      * \brief Returns the internal initializing stepper instance.
0411      * \return initializing_stepper
0412      */
0413
0414 } // odeint
0415 } // numeric
0416 } // boost
0417
0418
0419
0420 #endif // BOOST_NUMERIC_ODEINT_STEPPER_ADAMS_BASHFORTH_HPP_INCLUDED