stepper/base/symplectic_rkn_stepper_base.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp
0004
0005  [begin_description]
0006  Base class for symplectic Runge-Kutta-Nystrom steppers.
0007  [end_description]
0008
0009  Copyright 2011-2013 Karsten Ahnert
0010  Copyright 2011-2013 Mario Mulansky
0011  Copyright 2012 Christoph Koke
0012
0013  Distributed under the Boost Software License, Version 1.0.
0014  (See accompanying file LICENSE_1_0.txt or
0015  copy at http://www.boost.org/LICENSE_1_0.txt)
0016  */
0017
0018
0019 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_BASE_SYMPLECTIC_RKN_STEPPER_BASE_HPP_INCLUDED
0020 #define BOOST_NUMERIC_ODEINT_STEPPER_BASE_SYMPLECTIC_RKN_STEPPER_BASE_HPP_INCLUDED
0021
0022 #include <boost/array.hpp>
0023
0024 #include <boost/numeric/odeint/util/bind.hpp>
0025 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0026
0027 #include <boost/numeric/odeint/util/copy.hpp>
0028 #include <boost/numeric/odeint/util/is_pair.hpp>
0029
0030 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0031 #include <boost/numeric/odeint/util/resizer.hpp>
0032
0033 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0034
0035 #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
0036
0037
0038
0039
0040 namespace boost {
0041 namespace numeric {
0042 namespace odeint {
0043
0044
0045 template<
0046 size_t NumOfStages ,
0047 unsigned short Order ,
0048 class Coor ,
0049 class Momentum ,
0050 class Value ,
0051 class CoorDeriv ,
0052 class MomentumDeriv ,
0053 class Time ,
0054 class Algebra ,
0055 class Operations ,
0056 class Resizer
0057 >
0058 class symplectic_nystroem_stepper_base : public algebra_stepper_base< Algebra , Operations >
0059 {
0060
0061 public:
0062
0063     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
0064     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
0065     typedef typename algebra_stepper_base_type::operations_type operations_type;
0066
0067     const static size_t num_of_stages = NumOfStages;
0068     typedef Coor coor_type;
0069     typedef Momentum momentum_type;
0070     typedef std::pair< coor_type , momentum_type > state_type;
0071     typedef CoorDeriv coor_deriv_type;
0072     typedef state_wrapper< coor_deriv_type> wrapped_coor_deriv_type;
0073     typedef MomentumDeriv momentum_deriv_type;
0074     typedef state_wrapper< momentum_deriv_type > wrapped_momentum_deriv_type;
0075     typedef std::pair< coor_deriv_type , momentum_deriv_type > deriv_type;
0076     typedef Value value_type;
0077     typedef Time time_type;
0078     typedef Resizer resizer_type;
0079     typedef stepper_tag stepper_category;
0080
0081     #ifndef DOXYGEN_SKIP
0082     typedef symplectic_nystroem_stepper_base< NumOfStages , Order , Coor , Momentum , Value ,
0083             CoorDeriv , MomentumDeriv , Time , Algebra , Operations , Resizer > internal_stepper_base_type;
0084     #endif
0085     typedef unsigned short order_type;
0086
0087     static const order_type order_value = Order;
0088
0089     typedef boost::array< value_type , num_of_stages > coef_type;
0090
0091     symplectic_nystroem_stepper_base( const coef_type &coef_a , const coef_type &coef_b , const algebra_type &algebra = algebra_type() )
0092         : algebra_stepper_base_type( algebra ) , m_coef_a( coef_a ) , m_coef_b( coef_b ) ,
0093           m_dqdt_resizer() , m_dpdt_resizer() , m_dqdt() , m_dpdt()
0094     { }
0095
0096
0097     order_type order( void ) const
0098     {
0099         return order_value;
0100     }
0101
0102     /*
0103      * Version 1 : do_step( system , x , t , dt )
0104      *
0105      * This version does not solve the forwarding problem, boost.range can not be used.
0106      */
0107     template< class System , class StateInOut >
0108     void do_step( System system , const StateInOut &state , time_type t , time_type dt )
0109     {
0110         typedef typename odeint::unwrap_reference< System >::type system_type;
0111         do_step_impl( system , state , t , state , dt , typename is_pair< system_type >::type() );
0112     }
0113
0114     /**
0115      * \brief Same function as above. It differs only in a different const specifier in order
0116      * to solve the forwarding problem, can be used with Boost.Range.
0117      */
0118     template< class System , class StateInOut >
0119     void do_step( System system , StateInOut &state , time_type t , time_type dt )
0120     {
0121         typedef typename odeint::unwrap_reference< System >::type system_type;
0122         do_step_impl( system , state , t , state , dt , typename is_pair< system_type >::type() );
0123     }
0124
0125
0126
0127
0128     /*
0129      * Version 2 : do_step( system , q , p , t , dt );
0130      *
0131      * For Convenience
0132      *
0133      * The two overloads are needed in order to solve the forwarding problem.
0134      */
0135     template< class System , class CoorInOut , class MomentumInOut >
0136     void do_step( System system , CoorInOut &q , MomentumInOut &p , time_type t , time_type dt )
0137     {
0138         do_step( system , std::make_pair( detail::ref( q ) , detail::ref( p ) ) , t , dt );
0139     }
0140
0141     /**
0142      * \brief Same function as do_step( system , q , p , t , dt ). It differs only in a different const specifier in order
0143      * to solve the forwarding problem, can be called with Boost.Range.
0144      */
0145     template< class System , class CoorInOut , class MomentumInOut >
0146     void do_step( System system , const CoorInOut &q , const MomentumInOut &p , time_type t , time_type dt )
0147     {
0148         do_step( system , std::make_pair( detail::ref( q ) , detail::ref( p ) ) , t , dt );
0149     }
0150
0151
0152
0153
0154
0155     /*
0156      * Version 3 : do_step( system , in , t , out , dt )
0157      *
0158      * The forwarding problem is not solved in this version
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         typedef typename odeint::unwrap_reference< System >::type system_type;
0164         do_step_impl( system , in , t , out , dt , typename is_pair< system_type >::type() );
0165     }
0166
0167
0168     template< class StateType >
0169     void adjust_size( const StateType &x )
0170     {
0171         resize_dqdt( x );
0172         resize_dpdt( x );
0173     }
0174
0175     /** \brief Returns the coefficients a. */
0176     const coef_type& coef_a( void ) const { return m_coef_a; }
0177
0178     /** \brief Returns the coefficients b. */
0179     const coef_type& coef_b( void ) const { return m_coef_b; }
0180
0181 private:
0182
0183     // stepper for systems with function for dq/dt = f(p) and dp/dt = -f(q)
0184     template< class System , class StateIn , class StateOut >
0185     void do_step_impl( System system , const StateIn &in , time_type /* t */ , StateOut &out , time_type dt , boost::mpl::true_ )
0186     {
0187         typedef typename odeint::unwrap_reference< System >::type system_type;
0188         typedef typename odeint::unwrap_reference< typename system_type::first_type >::type coor_deriv_func_type;
0189         typedef typename odeint::unwrap_reference< typename system_type::second_type >::type momentum_deriv_func_type;
0190         system_type &sys = system;
0191         coor_deriv_func_type &coor_func = sys.first;
0192         momentum_deriv_func_type &momentum_func = sys.second;
0193
0194         typedef typename odeint::unwrap_reference< StateIn >::type state_in_type;
0195         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
0196         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
0197         const state_in_type &state_in = in;
0198         const coor_in_type &coor_in = state_in.first;
0199         const momentum_in_type &momentum_in = state_in.second;
0200
0201         typedef typename odeint::unwrap_reference< StateOut >::type state_out_type;
0202         typedef typename odeint::unwrap_reference< typename state_out_type::first_type >::type coor_out_type;
0203         typedef typename odeint::unwrap_reference< typename state_out_type::second_type >::type momentum_out_type;
0204         state_out_type &state_out = out;
0205         coor_out_type &coor_out = state_out.first;
0206         momentum_out_type &momentum_out = state_out.second;
0207
0208         m_dqdt_resizer.adjust_size( coor_in , detail::bind( &internal_stepper_base_type::template resize_dqdt< coor_in_type > , detail::ref( *this ) , detail::_1 ) );
0209         m_dpdt_resizer.adjust_size( momentum_in , detail::bind( &internal_stepper_base_type::template resize_dpdt< momentum_in_type > , detail::ref( *this ) , detail::_1 ) );
0210
0211         // ToDo: check sizes?
0212
0213         for( size_t l=0 ; l<num_of_stages ; ++l )
0214         {
0215             if( l == 0 )
0216             {
0217                 coor_func( momentum_in , m_dqdt.m_v );
0218                 this->m_algebra.for_each3( coor_out , coor_in , m_dqdt.m_v ,
0219                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
0220                 momentum_func( coor_out , m_dpdt.m_v );
0221                 this->m_algebra.for_each3( momentum_out , momentum_in , m_dpdt.m_v ,
0222                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
0223             }
0224             else
0225             {
0226                 coor_func( momentum_out , m_dqdt.m_v );
0227                 this->m_algebra.for_each3( coor_out , coor_out , m_dqdt.m_v ,
0228                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
0229                 momentum_func( coor_out , m_dpdt.m_v );
0230                 this->m_algebra.for_each3( momentum_out , momentum_out , m_dpdt.m_v ,
0231                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
0232             }
0233         }
0234     }
0235
0236
0237     // stepper for systems with only function dp /dt = -f(q), dq/dt = p, time not required but still expected for compatibility reasons
0238     template< class System , class StateIn , class StateOut >
0239     void do_step_impl( System system , const StateIn &in , time_type  /* t */ , StateOut &out , time_type dt , boost::mpl::false_ )
0240     {
0241         typedef typename odeint::unwrap_reference< System >::type momentum_deriv_func_type;
0242         momentum_deriv_func_type &momentum_func = system;
0243
0244         typedef typename odeint::unwrap_reference< StateIn >::type state_in_type;
0245         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
0246         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
0247         const state_in_type &state_in = in;
0248         const coor_in_type &coor_in = state_in.first;
0249         const momentum_in_type &momentum_in = state_in.second;
0250
0251         typedef typename odeint::unwrap_reference< StateOut >::type state_out_type;
0252         typedef typename odeint::unwrap_reference< typename state_out_type::first_type >::type coor_out_type;
0253         typedef typename odeint::unwrap_reference< typename state_out_type::second_type >::type momentum_out_type;
0254         state_out_type &state_out = out;
0255         coor_out_type &coor_out = state_out.first;
0256         momentum_out_type &momentum_out = state_out.second;
0257
0258
0259         // m_dqdt not required when called with momentum_func only - don't resize
0260         // m_dqdt_resizer.adjust_size( coor_in , detail::bind( &internal_stepper_base_type::template resize_dqdt< coor_in_type > , detail::ref( *this ) , detail::_1 ) );
0261         m_dpdt_resizer.adjust_size( momentum_in , detail::bind( &internal_stepper_base_type::template resize_dpdt< momentum_in_type > , detail::ref( *this ) , detail::_1 ) );
0262
0263
0264         // ToDo: check sizes?
0265
0266         // step 0
0267         this->m_algebra.for_each3( coor_out  , coor_in , momentum_in ,
0268                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[0] * dt ) );
0269         momentum_func( coor_out , m_dpdt.m_v );
0270         this->m_algebra.for_each3( momentum_out , momentum_in , m_dpdt.m_v ,
0271                                            typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[0] * dt ) );
0272
0273         for( size_t l=1 ; l<num_of_stages ; ++l )
0274         {
0275             this->m_algebra.for_each3( coor_out , coor_out , momentum_out ,
0276                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
0277             momentum_func( coor_out , m_dpdt.m_v );
0278             this->m_algebra.for_each3( momentum_out , momentum_out , m_dpdt.m_v ,
0279                                        typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
0280         }
0281     }
0282
0283     template< class StateIn >
0284     bool resize_dqdt( const StateIn &x )
0285     {
0286         return adjust_size_by_resizeability( m_dqdt , x , typename is_resizeable<coor_deriv_type>::type() );
0287     }
0288
0289     template< class StateIn >
0290     bool resize_dpdt( const StateIn &x )
0291     {
0292         return adjust_size_by_resizeability( m_dpdt , x , typename is_resizeable<momentum_deriv_type>::type() );
0293     }
0294
0295
0296     const coef_type m_coef_a;
0297     const coef_type m_coef_b;
0298
0299     resizer_type m_dqdt_resizer;
0300     resizer_type m_dpdt_resizer;
0301     wrapped_coor_deriv_type m_dqdt;
0302     wrapped_momentum_deriv_type m_dpdt;
0303
0304 };
0305
0306 /********* DOXYGEN *********/
0307
0308 /**
0309  * \class symplectic_nystroem_stepper_base
0310  * \brief Base class for all symplectic steppers of Nystroem type.
0311  *
0312  * This class is the base class for the symplectic Runge-Kutta-Nystroem steppers. Symplectic steppers are usually
0313  * used to solve Hamiltonian systems and they conserve the phase space volume, see
0314  * <a href="http://en.wikipedia.org/wiki/Symplectic_integrator">en.wikipedia.org/wiki/Symplectic_integrator</a>.
0315  * Furthermore, the energy is conserved
0316  * in average. In detail this class of steppers can be used to solve separable Hamiltonian systems which can be written
0317  * in the form H(q,p) = H1(p) + H2(q). q is usually called the coordinate, while p is the momentum. The equations of motion
0318  * are dq/dt = dH1/dp, dp/dt = -dH2/dq.
0319  *
0320  * ToDo : add formula for solver and explanation of the coefficients
0321  *
0322  * symplectic_nystroem_stepper_base uses odeints algebra and operation system. Step size and error estimation are not
0323  * provided for this class of solvers. It derives from algebra_stepper_base. Several `do_step` variants are provided:
0324  *
0325  * - `do_step( sys , x , t , dt )` - The classical `do_step` method. The sys can be either a pair of function objects
0326  *    for the coordinate or the momentum part or one function object for the momentum part. `x` is a pair of coordinate
0327  *    and momentum. The state is updated in-place.
0328  * - `do_step( sys , q , p , t , dt )` - This method is similar to the method above with the difference that the coordinate
0329  *    and the momentum are passed explicitly and not packed into a pair.
0330  * - `do_step( sys , x_in , t , x_out , dt )` - This method transforms the state out-of-place. `x_in` and `x_out` are here pairs
0331  *    of coordinate and momentum.
0332  *
0333  * \tparam NumOfStages Number of stages.
0334  * \tparam Order The order of the stepper.
0335  * \tparam Coor The type representing the coordinates q.
0336  * \tparam Momentum The type representing the coordinates p.
0337  * \tparam Value The basic value type. Should be something like float, double or a high-precision type.
0338  * \tparam CoorDeriv The type representing the time derivative of the coordinate dq/dt.
0339  * \tparam MomemtnumDeriv The type representing the time derivative of the momentum dp/dt.
0340  * \tparam Time The type representing the time t.
0341  * \tparam Algebra The algebra.
0342  * \tparam Operations The operations.
0343  * \tparam Resizer The resizer policy.
0344  */
0345
0346     /**
0347      * \fn symplectic_nystroem_stepper_base::symplectic_nystroem_stepper_base( const coef_type &coef_a , const coef_type &coef_b , const algebra_type &algebra )
0348      * \brief Constructs a symplectic_nystroem_stepper_base class. The parameters of the specific Nystroem method and the
0349      * algebra have to be passed.
0350      * \param coef_a The coefficients a.
0351      * \param coef_b The coefficients b.
0352      * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
0353      */
0354
0355     /**
0356      * \fn symplectic_nystroem_stepper_base::order( void ) const
0357      * \return Returns the order of the stepper.
0358      */
0359
0360     /**
0361      * \fn symplectic_nystroem_stepper_base::do_step( System system , const StateInOut &state , time_type t , time_type dt )
0362      * \brief This method performs one step. The system can be either a pair of two function object
0363      * describing the momentum part and the coordinate part or one function object describing only
0364      * the momentum part. In this case the coordinate is assumed to be trivial dq/dt = p. The state
0365      * is updated in-place.
0366      *
0367      * \note boost::ref or std::ref can be used for the system as well as for the state. So, it is correct
0368      * to write `stepper.do_step( make_pair( std::ref( fq ) , std::ref( fp ) ) , make_pair( std::ref( q ) , std::ref( p ) ) , t , dt )`.
0369      *
0370      * \note This method solves the forwarding problem.
0371      *
0372      * \param system The system, can be represented as a pair of two function object or one function object. See above.
0373      * \param state The state of the ODE. It is a pair of Coor and Momentum. The state is updated in-place, therefore, the
0374      * new value of the state will be written into this variable.
0375      * \param t The time of the ODE. It is not advanced by this method.
0376      * \param dt The time step.
0377      */
0378
0379     /**
0380      * \fn symplectic_nystroem_stepper_base::do_step( System system , CoorInOut &q , MomentumInOut &p , time_type t , time_type dt )
0381      * \brief This method performs one step. The system can be either a pair of two function object
0382      * describing the momentum part and the coordinate part or one function object describing only
0383      * the momentum part. In this case the coordinate is assumed to be trivial dq/dt = p. The state
0384      * is updated in-place.
0385      *
0386      * \note boost::ref or std::ref can be used for the system. So, it is correct
0387      * to write `stepper.do_step( make_pair( std::ref( fq ) , std::ref( fp ) ) , q , p , t , dt )`.
0388      *
0389      * \note This method solves the forwarding problem.
0390      *
0391      * \param system The system, can be represented as a pair of two function object or one function object. See above.
0392      * \param q The coordinate of the ODE. It is updated in-place. Therefore, the new value of the coordinate will be written
0393      * into this variable.
0394      * \param p The momentum of the ODE. It is updated in-place. Therefore, the new value of the momentum will be written info
0395      * this variable.
0396      * \param t The time of the ODE. It is not advanced by this method.
0397      * \param dt The time step.
0398      */
0399
0400     /**
0401      * \fn symplectic_nystroem_stepper_base::do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
0402      * \brief This method performs one step. The system can be either a pair of two function object
0403      * describing the momentum part and the coordinate part or one function object describing only
0404      * the momentum part. In this case the coordinate is assumed to be trivial dq/dt = p. The state
0405      * is updated out-of-place.
0406      *
0407      * \note boost::ref or std::ref can be used for the system. So, it is correct
0408      * to write `stepper.do_step( make_pair( std::ref( fq ) , std::ref( fp ) ) , x_in , t , x_out , dt )`.
0409      *
0410      * \note This method NOT solve the forwarding problem.
0411      *
0412      * \param system The system, can be represented as a pair of two function object or one function object. See above.
0413      * \param in The state of the ODE, which is a pair of coordinate and momentum. The state is updated out-of-place, therefore the
0414      * new value is written into out
0415      * \param t The time of the ODE. It is not advanced by this method.
0416      * \param out The new state of the ODE.
0417      * \param dt The time step.
0418      */
0419
0420     /**
0421      * \fn symplectic_nystroem_stepper_base::adjust_size( const StateType &x )
0422      * \brief Adjust the size of all temporaries in the stepper manually.
0423      * \param x A state from which the size of the temporaries to be resized is deduced.
0424      */
0425
0426 } // namespace odeint
0427 } // namespace numeric
0428 } // namespace boost
0429
0430
0431 #endif // BOOST_NUMERIC_ODEINT_STEPPER_BASE_SYMPLECTIC_RKN_STEPPER_BASE_HPP_INCLUDED