odeint/stepper/modified_midpoint.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/stepper/modified_midpoint.hpp
0004
0005  [begin_description]
0006  Modified midpoint method for the use in Burlish-Stoer stepper.
0007  [end_description]
0008
0009  Copyright 2011-2013 Mario Mulansky
0010  Copyright 2011-2013 Karsten Ahnert
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_MODIFIED_MIDPOINT_HPP_INCLUDED
0020 #define BOOST_NUMERIC_ODEINT_STEPPER_MODIFIED_MIDPOINT_HPP_INCLUDED
0021
0022 #include <vector>
0023
0024 #include <boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp>
0025 #include <boost/numeric/odeint/util/resizer.hpp>
0026 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0027 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
0028 #include <boost/numeric/odeint/algebra/default_operations.hpp>
0029 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
0030 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
0031 #include <boost/numeric/odeint/util/copy.hpp>
0032
0033 namespace boost {
0034 namespace numeric {
0035 namespace odeint {
0036
0037 template<
0038 class State ,
0039 class Value = double ,
0040 class Deriv = State ,
0041 class Time = Value ,
0042 class Algebra = typename algebra_dispatcher< State >::algebra_type ,
0043 class Operations = typename operations_dispatcher< State >::operations_type ,
0044 class Resizer = initially_resizer
0045 >
0046 #ifndef DOXYGEN_SKIP
0047 class modified_midpoint
0048 : public explicit_stepper_base<
0049   modified_midpoint< State , Value , Deriv , Time , Algebra , Operations , Resizer > ,
0050   2 , State , Value , Deriv , Time , Algebra , Operations , Resizer >
0051 #else
0052 class modified_midpoint : public explicit_stepper_base
0053 #endif
0054 {
0055
0056 public :
0057
0058     typedef explicit_stepper_base<
0059     modified_midpoint< State , Value , Deriv , Time , Algebra , Operations , Resizer > ,
0060     2 , State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_base_type;
0061
0062     typedef typename stepper_base_type::state_type state_type;
0063     typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
0064     typedef typename stepper_base_type::value_type value_type;
0065     typedef typename stepper_base_type::deriv_type deriv_type;
0066     typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
0067     typedef typename stepper_base_type::time_type time_type;
0068     typedef typename stepper_base_type::algebra_type algebra_type;
0069     typedef typename stepper_base_type::operations_type operations_type;
0070     typedef typename stepper_base_type::resizer_type resizer_type;
0071     typedef typename stepper_base_type::stepper_type stepper_type;
0072
0073
0074     modified_midpoint( unsigned short steps = 2 , const algebra_type &algebra = algebra_type() )
0075     : stepper_base_type( algebra ) , m_steps( steps )
0076     { }
0077
0078     template< class System , class StateIn , class DerivIn , class StateOut >
0079     void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
0080     {
0081         static const value_type val1 = static_cast< value_type >( 1 );
0082         static const value_type val05 = static_cast< value_type >( 1 ) / static_cast< value_type >( 2 );
0083
0084         m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
0085
0086         const time_type h = dt /  static_cast<value_type>( m_steps );
0087         const time_type h2 = static_cast<value_type>(2) * h;
0088
0089         typename odeint::unwrap_reference< System >::type &sys = system;
0090
0091         time_type th = t + h;
0092
0093         // m_x1 = x + h*dxdt
0094         stepper_base_type::m_algebra.for_each3( m_x1.m_v , in , dxdt ,
0095                 typename operations_type::template scale_sum2< value_type , time_type >( val1 , h ) );
0096
0097         sys( m_x1.m_v , m_dxdt.m_v , th );
0098
0099         boost::numeric::odeint::copy( in , m_x0.m_v );
0100
0101         unsigned short i = 1;
0102         while( i != m_steps )
0103         {
0104             // general step
0105             //tmp = m_x1; m_x1 = m_x0 + h2*m_dxdt; m_x0 = tmp
0106             stepper_base_type::m_algebra.for_each3( m_x1.m_v , m_x0.m_v , m_dxdt.m_v ,
0107                     typename operations_type::template scale_sum_swap2< value_type , time_type >( val1 , h2 ) );
0108             th += h;
0109             sys( m_x1.m_v , m_dxdt.m_v , th);
0110             i++;
0111         }
0112
0113         // last step
0114         // x = 0.5*( m_x0 + m_x1 + h*m_dxdt )
0115         stepper_base_type::m_algebra.for_each4( out , m_x0.m_v , m_x1.m_v , m_dxdt.m_v ,
0116                 typename operations_type::template scale_sum3< value_type , value_type , time_type >( val05 , val05 , val05*h ) );
0117     }
0118
0119
0120     void set_steps( unsigned short steps )
0121     {   m_steps = steps; }
0122
0123
0124     unsigned short steps( void ) const
0125     {   return m_steps; }
0126
0127
0128     template< class StateIn >
0129     void adjust_size( const StateIn &x )
0130     {
0131         resize_impl( x );
0132         stepper_base_type::adjust_size( x );
0133     }
0134
0135 private:
0136
0137     template< class StateIn >
0138     bool resize_impl( const StateIn &x )
0139     {
0140         bool resized( false );
0141         resized |= adjust_size_by_resizeability( m_x0 , x , typename is_resizeable<state_type>::type() );
0142         resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() );
0143         resized |= adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
0144         return resized;
0145     }
0146
0147
0148     unsigned short m_steps;
0149
0150     resizer_type m_resizer;
0151
0152     wrapped_state_type m_x0;
0153     wrapped_state_type m_x1;
0154     wrapped_deriv_type m_dxdt;
0155
0156 };
0157
0158
0159 /* Modified midpoint which stores derivatives and state at dt/2 in some external storage for later usage in dense output calculation
0160  * This Stepper is for use in Bulirsch Stoer only. It DOES NOT meet any stepper concept.
0161  */
0162 template<
0163 class State ,
0164 class Value = double ,
0165 class Deriv = State ,
0166 class Time = Value ,
0167 class Algebra = typename algebra_dispatcher< State >::algebra_type ,
0168 class Operations = typename operations_dispatcher< State >::operations_type ,
0169 class Resizer = initially_resizer
0170 >
0171 class modified_midpoint_dense_out
0172 {
0173
0174 public :
0175
0176     typedef State state_type;
0177     typedef Value value_type;
0178     typedef Deriv deriv_type;
0179     typedef Time time_type;
0180     typedef Algebra algebra_type;
0181     typedef Operations operations_type;
0182     typedef Resizer resizer_type;
0183     typedef state_wrapper< state_type > wrapped_state_type;
0184     typedef state_wrapper< deriv_type > wrapped_deriv_type;
0185
0186     typedef modified_midpoint_dense_out< State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_type;
0187     typedef std::vector< wrapped_deriv_type > deriv_table_type;
0188
0189     modified_midpoint_dense_out( unsigned short steps = 2 , const algebra_type &algebra = algebra_type() )
0190     : m_algebra( algebra ) , m_steps( steps )
0191     { }
0192
0193     /*
0194      * performs a modified midpoint step with m_steps intermediate points
0195      * stores approximation for x(t+dt/2) in x_mp and all evaluated function results in derivs
0196      *
0197      */
0198
0199     template< class System , class StateIn , class DerivIn , class StateOut >
0200     void do_step( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt ,
0201             state_type &x_mp , deriv_table_type &derivs )
0202     {
0203
0204         static const value_type val1 = static_cast< value_type >( 1 );
0205         static const value_type val05 = static_cast< value_type >( 1 ) / static_cast< value_type >( 2 );
0206
0207         m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize< StateIn > , detail::ref( *this ) , detail::_1 ) );
0208
0209         const time_type h = dt / static_cast<value_type>( m_steps );
0210         const time_type h2 = static_cast<value_type>( 2 ) * h;
0211
0212         typename odeint::unwrap_reference< System >::type &sys = system;
0213
0214         time_type th = t + h;
0215
0216         // m_x1 = x + h*dxdt
0217         m_algebra.for_each3( m_x1.m_v , in , dxdt ,
0218                 typename operations_type::template scale_sum2< value_type , time_type >( val1 , h ) );
0219
0220         if( m_steps == 2 )
0221             // result of first step already gives approximation at the center of the interval
0222             boost::numeric::odeint::copy( m_x1.m_v , x_mp );
0223
0224         sys( m_x1.m_v , derivs[0].m_v , th );
0225
0226         boost::numeric::odeint::copy( in , m_x0.m_v );
0227
0228         unsigned short i = 1;
0229         while( i != m_steps )
0230         {
0231             // general step
0232             //tmp = m_x1; m_x1 = m_x0 + h2*m_dxdt; m_x0 = tmp
0233             m_algebra.for_each3( m_x1.m_v , m_x0.m_v , derivs[i-1].m_v ,
0234                     typename operations_type::template scale_sum_swap2< value_type , time_type >( val1 , h2 ) );
0235             if( i == m_steps/2-1 )
0236                 // save approximation at the center of the interval
0237                 boost::numeric::odeint::copy( m_x1.m_v , x_mp );
0238
0239             th += h;
0240             sys( m_x1.m_v , derivs[i].m_v , th);
0241             i++;
0242         }
0243
0244         // last step
0245         // x = 0.5*( m_x0 + m_x1 + h*m_dxdt )
0246         m_algebra.for_each4( out , m_x0.m_v , m_x1.m_v , derivs[m_steps-1].m_v ,
0247                 typename operations_type::template scale_sum3< value_type , value_type , time_type >( val05 , val05 , val05*h ) );
0248     }
0249
0250
0251     void set_steps( unsigned short steps )
0252     {   m_steps = steps; }
0253
0254
0255     unsigned short steps( void ) const
0256     {   return m_steps; }
0257
0258
0259     template< class StateIn >
0260     bool resize( const StateIn &x )
0261     {
0262         bool resized( false );
0263         resized |= adjust_size_by_resizeability( m_x0 , x , typename is_resizeable<state_type>::type() );
0264         resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() );
0265         return resized;
0266     }
0267
0268     template< class StateIn >
0269     void adjust_size( const StateIn &x )
0270     {
0271         resize( x );
0272     }
0273
0274 private:
0275
0276     algebra_type m_algebra;
0277
0278     unsigned short m_steps;
0279
0280     resizer_type m_resizer;
0281
0282     wrapped_state_type m_x0;
0283     wrapped_state_type m_x1;
0284
0285 };
0286
0287
0288
0289 /********** DOXYGEN ***********/
0290
0291 /**
0292  * \class modified_midpoint
0293  *
0294  * Implementation of the modified midpoint method with a configurable
0295  * number of intermediate steps. This class is used by the Bulirsch-Stoer
0296  * algorithm and is not meant for direct usage.
0297  */
0298
0299
0300 /**
0301  * \class modified_midpoint_dense_out
0302  *
0303  * Implementation of the modified midpoint method with a configurable
0304  * number of intermediate steps. This class is used by the dense output
0305  * Bulirsch-Stoer algorithm and is not meant for direct usage.
0306  * \note This stepper is for internal use only and does not meet
0307  * any stepper concept.
0308  */
0309
0310
0311 }
0312 }
0313 }
0314
0315 #endif // BOOST_NUMERIC_ODEINT_STEPPER_MODIFIED_MIDPOINT_HPP_INCLUDED