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 <array>
0023 #include <type_traits>
0024
0025 #include <boost/numeric/odeint/util/bind.hpp>
0026 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0027
0028 #include <boost/numeric/odeint/util/copy.hpp>
0029 #include <boost/numeric/odeint/util/is_pair.hpp>
0030
0031 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0032 #include <boost/numeric/odeint/util/resizer.hpp>
0033
0034 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0035
0036 #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
0037
0038
0039
0040
0041 namespace boost {
0042 namespace numeric {
0043 namespace odeint {
0044
0045
0046 template<
0047 size_t NumOfStages ,
0048 unsigned short Order ,
0049 class Coor ,
0050 class Momentum ,
0051 class Value ,
0052 class CoorDeriv ,
0053 class MomentumDeriv ,
0054 class Time ,
0055 class Algebra ,
0056 class Operations ,
0057 class Resizer
0058 >
0059 class symplectic_nystroem_stepper_base : public algebra_stepper_base< Algebra , Operations >
0060 {
0061
0062 public:
0063
0064     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
0065     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
0066     typedef typename algebra_stepper_base_type::operations_type operations_type;
0067
0068     const static size_t num_of_stages = NumOfStages;
0069     typedef Coor coor_type;
0070     typedef Momentum momentum_type;
0071     typedef std::pair< coor_type , momentum_type > state_type;
0072     typedef CoorDeriv coor_deriv_type;
0073     typedef state_wrapper< coor_deriv_type> wrapped_coor_deriv_type;
0074     typedef MomentumDeriv momentum_deriv_type;
0075     typedef state_wrapper< momentum_deriv_type > wrapped_momentum_deriv_type;
0076     typedef std::pair< coor_deriv_type , momentum_deriv_type > deriv_type;
0077     typedef Value value_type;
0078     typedef Time time_type;
0079     typedef Resizer resizer_type;
0080     typedef stepper_tag stepper_category;
0081
0082     #ifndef DOXYGEN_SKIP
0083     typedef symplectic_nystroem_stepper_base< NumOfStages , Order , Coor , Momentum , Value ,
0084             CoorDeriv , MomentumDeriv , Time , Algebra , Operations , Resizer > internal_stepper_base_type;
0085     #endif
0086     typedef unsigned short order_type;
0087
0088     static const order_type order_value = Order;
0089
0090     typedef std::array< value_type , num_of_stages > coef_type;
0091
0092     symplectic_nystroem_stepper_base( const coef_type &coef_a , const coef_type &coef_b , const algebra_type &algebra = algebra_type() )
0093         : algebra_stepper_base_type( algebra ) , m_coef_a( coef_a ) , m_coef_b( coef_b ) ,
0094           m_dqdt_resizer() , m_dpdt_resizer() , m_dqdt() , m_dpdt()
0095     { }
0096
0097
0098     order_type order( void ) const
0099     {
0100         return order_value;
0101     }
0102
0103     /*
0104      * Version 1 : do_step( system , x , t , dt )
0105      *
0106      * This version does not solve the forwarding problem, boost.range can not be used.
0107      */
0108     template< class System , class StateInOut >
0109     void do_step( System system , const StateInOut &state , time_type t , time_type dt )
0110     {
0111         typedef typename odeint::unwrap_reference< System >::type system_type;
0112         do_step_impl( system , state , t , state , dt , typename is_pair< system_type >::type() );
0113     }
0114
0115     /**
0116      * \brief Same function as above. It differs only in a different const specifier in order
0117      * to solve the forwarding problem, can be used with Boost.Range.
0118      */
0119     template< class System , class StateInOut >
0120     void do_step( System system , StateInOut &state , time_type t , time_type dt )
0121     {
0122         typedef typename odeint::unwrap_reference< System >::type system_type;
0123         do_step_impl( system , state , t , state , dt , typename is_pair< system_type >::type() );
0124     }
0125
0126
0127
0128
0129     /*
0130      * Version 2 : do_step( system , q , p , t , dt );
0131      *
0132      * For Convenience
0133      *
0134      * The two overloads are needed in order to solve the forwarding problem.
0135      */
0136     template< class System , class CoorInOut , class MomentumInOut >
0137     void do_step( System system , CoorInOut &q , MomentumInOut &p , time_type t , time_type dt )
0138     {
0139         do_step( system , std::make_pair( std::ref( q ) , std::ref( p ) ) , t , dt );
0140     }
0141
0142     /**
0143      * \brief Same function as do_step( system , q , p , t , dt ). It differs only in a different const specifier in order
0144      * to solve the forwarding problem, can be called with Boost.Range.
0145      */
0146     template< class System , class CoorInOut , class MomentumInOut >
0147     void do_step( System system , const CoorInOut &q , const MomentumInOut &p , time_type t , time_type dt )
0148     {
0149         do_step( system , std::make_pair( std::ref( q ) , std::ref( p ) ) , t , dt );
0150     }
0151
0152
0153
0154
0155
0156     /*
0157      * Version 3 : do_step( system , in , t , out , dt )
0158      *
0159      * The forwarding problem is not solved in this version
0160      */
0161     template< class System , class StateIn , class StateOut >
0162     void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
0163     {
0164         typedef typename odeint::unwrap_reference< System >::type system_type;
0165         do_step_impl( system , in , t , out , dt , typename is_pair< system_type >::type() );
0166     }
0167
0168
0169     template< class StateType >
0170     void adjust_size( const StateType &x )
0171     {
0172         resize_dqdt( x );
0173         resize_dpdt( x );
0174     }
0175
0176     /** \brief Returns the coefficients a. */
0177     const coef_type& coef_a( void ) const { return m_coef_a; }
0178
0179     /** \brief Returns the coefficients b. */
0180     const coef_type& coef_b( void ) const { return m_coef_b; }
0181
0182 private:
0183
0184     // stepper for systems with function for dq/dt = f(p) and dp/dt = -f(q)
0185     template< class System , class StateIn , class StateOut >
0186     void do_step_impl( System system , const StateIn &in , time_type /* t */ , StateOut &out , time_type dt , std::integral_constant<bool, true> )
0187     {
0188         typedef typename odeint::unwrap_reference< System >::type system_type;
0189         typedef typename odeint::unwrap_reference< typename system_type::first_type >::type coor_deriv_func_type;
0190         typedef typename odeint::unwrap_reference< typename system_type::second_type >::type momentum_deriv_func_type;
0191         system_type &sys = system;
0192         coor_deriv_func_type &coor_func = sys.first;
0193         momentum_deriv_func_type &momentum_func = sys.second;
0194
0195         typedef typename odeint::unwrap_reference< StateIn >::type state_in_type;
0196         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
0197         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
0198         const state_in_type &state_in = in;
0199         const coor_in_type &coor_in = state_in.first;
0200         const momentum_in_type &momentum_in = state_in.second;
0201
0202         typedef typename odeint::unwrap_reference< StateOut >::type state_out_type;
0203         typedef typename odeint::unwrap_reference< typename state_out_type::first_type >::type coor_out_type;
0204         typedef typename odeint::unwrap_reference< typename state_out_type::second_type >::type momentum_out_type;
0205         state_out_type &state_out = out;
0206         coor_out_type &coor_out = state_out.first;
0207         momentum_out_type &momentum_out = state_out.second;
0208
0209         m_dqdt_resizer.adjust_size(coor_in, [this](auto&& arg) { return this->resize_dqdt<coor_in_type>(std::forward<decltype(arg)>(arg)); });
0210         m_dpdt_resizer.adjust_size(momentum_in, [this](auto&& arg) { return this->resize_dpdt<momentum_in_type>(std::forward<decltype(arg)>(arg)); });
0211
0212         // ToDo: check sizes?
0213
0214         for( size_t l=0 ; l<num_of_stages ; ++l )
0215         {
0216             if( l == 0 )
0217             {
0218                 coor_func( momentum_in , m_dqdt.m_v );
0219                 this->m_algebra.for_each3( coor_out , coor_in , m_dqdt.m_v ,
0220                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
0221                 momentum_func( coor_out , m_dpdt.m_v );
0222                 this->m_algebra.for_each3( momentum_out , momentum_in , m_dpdt.m_v ,
0223                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
0224             }
0225             else
0226             {
0227                 coor_func( momentum_out , m_dqdt.m_v );
0228                 this->m_algebra.for_each3( coor_out , coor_out , m_dqdt.m_v ,
0229                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
0230                 momentum_func( coor_out , m_dpdt.m_v );
0231                 this->m_algebra.for_each3( momentum_out , momentum_out , m_dpdt.m_v ,
0232                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
0233             }
0234         }
0235     }
0236
0237
0238     // stepper for systems with only function dp /dt = -f(q), dq/dt = p, time not required but still expected for compatibility reasons
0239     template< class System , class StateIn , class StateOut >
0240     void do_step_impl( System system , const StateIn &in , time_type  /* t */ , StateOut &out , time_type dt , std::integral_constant<bool, false> )
0241     {
0242         typedef typename odeint::unwrap_reference< System >::type momentum_deriv_func_type;
0243         momentum_deriv_func_type &momentum_func = system;
0244
0245         typedef typename odeint::unwrap_reference< StateIn >::type state_in_type;
0246         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
0247         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
0248         const state_in_type &state_in = in;
0249         const coor_in_type &coor_in = state_in.first;
0250         const momentum_in_type &momentum_in = state_in.second;
0251
0252         typedef typename odeint::unwrap_reference< StateOut >::type state_out_type;
0253         typedef typename odeint::unwrap_reference< typename state_out_type::first_type >::type coor_out_type;
0254         typedef typename odeint::unwrap_reference< typename state_out_type::second_type >::type momentum_out_type;
0255         state_out_type &state_out = out;
0256         coor_out_type &coor_out = state_out.first;
0257         momentum_out_type &momentum_out = state_out.second;
0258
0259
0260         // m_dqdt not required when called with momentum_func only - don't resize
0261         m_dpdt_resizer.adjust_size(momentum_in, [this](auto&& arg) { return this->resize_dpdt<momentum_in_type>(std::forward<decltype(arg)>(arg)); });
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