odeint/stepper/dense_output_runge_kutta.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp
0004
0005  [begin_description]
0006  Implementation of the Dense-output stepper for all steppers. Note, that this class does
0007  not computes the result but serves as an interface.
0008  [end_description]
0009
0010  Copyright 2011-2013 Karsten Ahnert
0011  Copyright 2011-2015 Mario Mulansky
0012  Copyright 2012 Christoph Koke
0013
0014  Distributed under the Boost Software License, Version 1.0.
0015  (See accompanying file LICENSE_1_0.txt or
0016  copy at http://www.boost.org/LICENSE_1_0.txt)
0017  */
0018
0019
0020 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED
0021 #define BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED
0022
0023
0024 #include <utility>
0025 #include <stdexcept>
0026
0027 #include <boost/throw_exception.hpp>
0028
0029 #include <boost/numeric/odeint/util/bind.hpp>
0030
0031 #include <boost/numeric/odeint/util/copy.hpp>
0032
0033 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0034 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0035 #include <boost/numeric/odeint/util/resizer.hpp>
0036
0037 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
0038 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0039
0040 #include <boost/numeric/odeint/integrate/max_step_checker.hpp>
0041
0042 namespace boost {
0043 namespace numeric {
0044 namespace odeint {
0045
0046 template< class Stepper , class StepperCategory = typename Stepper::stepper_category >
0047 class dense_output_runge_kutta;
0048
0049
0050 /**
0051  * \brief The class representing dense-output Runge-Kutta steppers.
0052  * \note In this stepper, the initialize method has to be called before using
0053  * the do_step method.
0054  *
0055  * The dense-output functionality allows to interpolate the solution between
0056  * subsequent integration points using intermediate results obtained during the
0057  * computation. This version works based on a normal stepper without step-size
0058  * control.
0059  *
0060  *
0061  * \tparam Stepper The stepper type of the underlying algorithm.
0062  */
0063 template< class Stepper >
0064 class dense_output_runge_kutta< Stepper , stepper_tag >
0065 {
0066
0067 public:
0068
0069     /*
0070      * We do not need all typedefs.
0071      */
0072     typedef Stepper stepper_type;
0073     typedef typename stepper_type::state_type state_type;
0074     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
0075     typedef typename stepper_type::value_type value_type;
0076     typedef typename stepper_type::deriv_type deriv_type;
0077     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
0078     typedef typename stepper_type::time_type time_type;
0079     typedef typename stepper_type::algebra_type algebra_type;
0080     typedef typename stepper_type::operations_type operations_type;
0081     typedef typename stepper_type::resizer_type resizer_type;
0082     typedef dense_output_stepper_tag stepper_category;
0083     typedef dense_output_runge_kutta< Stepper > dense_output_stepper_type;
0084
0085
0086     /**
0087      * \brief Constructs the dense_output_runge_kutta class. An instance of the
0088      * underlying stepper can be provided.
0089      * \param stepper An instance of the underlying stepper.
0090      */
0091     dense_output_runge_kutta( const stepper_type &stepper = stepper_type() )
0092     : m_stepper( stepper ) , m_resizer() ,
0093       m_x1() , m_x2() , m_current_state_x1( true ) ,
0094       m_t() , m_t_old() , m_dt()
0095     { }
0096
0097
0098     /**
0099      * \brief Initializes the stepper. Has to be called before do_step can be
0100      * used to set the initial conditions and the step size.
0101      * \param x0 The initial state of the ODE which should be solved.
0102      * \param t0 The initial time, at which the step should be performed.
0103      * \param dt0 The step size.
0104      */
0105     template< class StateType >
0106     void initialize( const StateType &x0 , time_type t0 , time_type dt0 )
0107     {
0108         m_resizer.adjust_size( x0 , detail::bind( &dense_output_stepper_type::template resize_impl< StateType > , detail::ref( *this ) , detail::_1 ) );
0109         boost::numeric::odeint::copy( x0 , get_current_state() );
0110         m_t = t0;
0111         m_dt = dt0;
0112     }
0113
0114     /**
0115      * \brief Does one time step.
0116      * \note initialize has to be called before using this method to set the
0117      * initial conditions x,t and the stepsize.
0118      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
0119      *               Simple System concept.
0120      * \return Pair with start and end time of the integration step.
0121      */
0122     template< class System >
0123     std::pair< time_type , time_type > do_step( System system )
0124     {
0125         m_stepper.do_step( system , get_current_state() , m_t , get_old_state() , m_dt );
0126         m_t_old = m_t;
0127         m_t += m_dt;
0128         toggle_current_state();
0129         return std::make_pair( m_t_old , m_t );
0130     }
0131
0132     /*
0133      * The next two overloads are needed to solve the forwarding problem
0134      */
0135
0136     /**
0137      * \brief Calculates the solution at an intermediate point.
0138      * \param t The time at which the solution should be calculated, has to be
0139      * in the current time interval.
0140      * \param x The output variable where the result is written into.
0141      */
0142     template< class StateOut >
0143     void calc_state( time_type t , StateOut &x ) const
0144     {
0145         if( t == current_time() )
0146         {
0147             boost::numeric::odeint::copy( get_current_state() , x );
0148         }
0149         m_stepper.calc_state( x , t , get_old_state() , m_t_old , get_current_state() , m_t );
0150     }
0151
0152     /**
0153      * \brief Calculates the solution at an intermediate point. Solves the forwarding problem
0154      * \param t The time at which the solution should be calculated, has to be
0155      * in the current time interval.
0156      * \param x The output variable where the result is written into, can be a boost range.
0157      */
0158     template< class StateOut >
0159     void calc_state( time_type t , const StateOut &x ) const
0160     {
0161         m_stepper.calc_state( x , t , get_old_state() , m_t_old , get_current_state() , m_t );
0162     }
0163
0164     /**
0165      * \brief Adjust the size of all temporaries in the stepper manually.
0166      * \param x A state from which the size of the temporaries to be resized is deduced.
0167      */
0168     template< class StateType >
0169     void adjust_size( const StateType &x )
0170     {
0171         resize_impl( x );
0172         m_stepper.stepper().resize( x );
0173     }
0174
0175     /**
0176      * \brief Returns the current state of the solution.
0177      * \return The current state of the solution x(t).
0178      */
0179     const state_type& current_state( void ) const
0180     {
0181         return get_current_state();
0182     }
0183
0184     /**
0185      * \brief Returns the current time of the solution.
0186      * \return The current time of the solution t.
0187      */
0188     time_type current_time( void ) const
0189     {
0190         return m_t;
0191     }
0192
0193     /**
0194      * \brief Returns the last state of the solution.
0195      * \return The last state of the solution x(t-dt).
0196      */
0197     const state_type& previous_state( void ) const
0198     {
0199         return get_old_state();
0200     }
0201
0202     /**
0203      * \brief Returns the last time of the solution.
0204      * \return The last time of the solution t-dt.
0205      */
0206     time_type previous_time( void ) const
0207     {
0208         return m_t_old;
0209     }
0210
0211     /**
0212      * \brief Returns the current time step.
0213      * \return dt.
0214      */
0215     time_type current_time_step( void ) const
0216     {
0217         return m_dt;
0218     }
0219
0220
0221 private:
0222
0223     state_type& get_current_state( void )
0224     {
0225         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
0226     }
0227
0228     const state_type& get_current_state( void ) const
0229     {
0230         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
0231     }
0232
0233     state_type& get_old_state( void )
0234     {
0235         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
0236     }
0237
0238     const state_type& get_old_state( void ) const
0239     {
0240         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
0241     }
0242
0243     void toggle_current_state( void )
0244     {
0245         m_current_state_x1 = ! m_current_state_x1;
0246     }
0247
0248
0249     template< class StateIn >
0250     bool resize_impl( const StateIn &x )
0251     {
0252         bool resized = false;
0253         resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() );
0254         resized |= adjust_size_by_resizeability( m_x2 , x , typename is_resizeable<state_type>::type() );
0255         return resized;
0256     }
0257
0258
0259     stepper_type m_stepper;
0260     resizer_type m_resizer;
0261     wrapped_state_type m_x1 , m_x2;
0262     bool m_current_state_x1;    // if true, the current state is m_x1
0263     time_type m_t , m_t_old , m_dt;
0264
0265 };
0266
0267
0268
0269
0270
0271 /**
0272  * \brief The class representing dense-output Runge-Kutta steppers with FSAL property.
0273  *
0274  * The interface is the same as for dense_output_runge_kutta< Stepper , stepper_tag >.
0275  * This class provides dense output functionality based on methods with step size controlled
0276  *
0277  *
0278  * \tparam Stepper The stepper type of the underlying algorithm.
0279  */
0280 template< class Stepper >
0281 class dense_output_runge_kutta< Stepper , explicit_controlled_stepper_fsal_tag >
0282 {
0283 public:
0284
0285     /*
0286      * We do not need all typedefs.
0287      */
0288     typedef Stepper controlled_stepper_type;
0289
0290     typedef typename controlled_stepper_type::stepper_type stepper_type;
0291     typedef typename stepper_type::state_type state_type;
0292     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
0293     typedef typename stepper_type::value_type value_type;
0294     typedef typename stepper_type::deriv_type deriv_type;
0295     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
0296     typedef typename stepper_type::time_type time_type;
0297     typedef typename stepper_type::algebra_type algebra_type;
0298     typedef typename stepper_type::operations_type operations_type;
0299     typedef typename stepper_type::resizer_type resizer_type;
0300     typedef dense_output_stepper_tag stepper_category;
0301     typedef dense_output_runge_kutta< Stepper > dense_output_stepper_type;
0302
0303
0304     dense_output_runge_kutta( const controlled_stepper_type &stepper = controlled_stepper_type() )
0305     : m_stepper( stepper ) , m_resizer() ,
0306       m_current_state_x1( true ) ,
0307       m_x1() , m_x2() , m_dxdt1() , m_dxdt2() ,
0308       m_t() , m_t_old() , m_dt() ,
0309       m_is_deriv_initialized( false )
0310     { }
0311
0312
0313     template< class StateType >
0314     void initialize( const StateType &x0 , time_type t0 , time_type dt0 )
0315     {
0316         m_resizer.adjust_size( x0 , detail::bind( &dense_output_stepper_type::template resize< StateType > , detail::ref( *this ) , detail::_1 ) );
0317         boost::numeric::odeint::copy( x0 , get_current_state() );
0318         m_t = t0;
0319         m_dt = dt0;
0320         m_is_deriv_initialized = false;
0321     }
0322
0323     template< class System >
0324     std::pair< time_type , time_type > do_step( System system )
0325     {
0326         if( !m_is_deriv_initialized )
0327         {
0328             typename odeint::unwrap_reference< System >::type &sys = system;
0329             sys( get_current_state() , get_current_deriv() , m_t );
0330             m_is_deriv_initialized = true;
0331         }
0332
0333         failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
0334         controlled_step_result res = fail;
0335         m_t_old = m_t;
0336         do
0337         {
0338             res = m_stepper.try_step( system , get_current_state() , get_current_deriv() , m_t ,
0339                                       get_old_state() , get_old_deriv() , m_dt );
0340             fail_checker();  // check for overflow of failed steps
0341         }
0342         while( res == fail );
0343         toggle_current_state();
0344         return std::make_pair( m_t_old , m_t );
0345     }
0346
0347
0348     /*
0349      * The two overloads are needed in order to solve the forwarding problem.
0350      */
0351     template< class StateOut >
0352     void calc_state( time_type t , StateOut &x ) const
0353     {
0354         m_stepper.stepper().calc_state( t , x , get_old_state() , get_old_deriv() , m_t_old ,
0355                                         get_current_state() , get_current_deriv() , m_t );
0356     }
0357
0358     template< class StateOut >
0359     void calc_state( time_type t , const StateOut &x ) const
0360     {
0361         m_stepper.stepper().calc_state( t , x , get_old_state() , get_old_deriv() , m_t_old ,
0362                                         get_current_state() , get_current_deriv() , m_t );
0363     }
0364
0365
0366     template< class StateIn >
0367     bool resize( const StateIn &x )
0368     {
0369         bool resized = false;
0370         resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() );
0371         resized |= adjust_size_by_resizeability( m_x2 , x , typename is_resizeable<state_type>::type() );
0372         resized |= adjust_size_by_resizeability( m_dxdt1 , x , typename is_resizeable<deriv_type>::type() );
0373         resized |= adjust_size_by_resizeability( m_dxdt2 , x , typename is_resizeable<deriv_type>::type() );
0374         return resized;
0375     }
0376
0377
0378     template< class StateType >
0379     void adjust_size( const StateType &x )
0380     {
0381         resize( x );
0382         m_stepper.stepper().resize( x );
0383     }
0384
0385     const state_type& current_state( void ) const
0386     {
0387         return get_current_state();
0388     }
0389
0390     time_type current_time( void ) const
0391     {
0392         return m_t;
0393     }
0394
0395     const state_type& previous_state( void ) const
0396     {
0397         return get_old_state();
0398     }
0399
0400     time_type previous_time( void ) const
0401     {
0402         return m_t_old;
0403     }
0404
0405     time_type current_time_step( void ) const
0406     {
0407         return m_dt;
0408     }
0409
0410
0411 private:
0412
0413     state_type& get_current_state( void )
0414     {
0415         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
0416     }
0417
0418     const state_type& get_current_state( void ) const
0419     {
0420         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
0421     }
0422
0423     state_type& get_old_state( void )
0424     {
0425         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
0426     }
0427
0428     const state_type& get_old_state( void ) const
0429     {
0430         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
0431     }
0432
0433     deriv_type& get_current_deriv( void )
0434     {
0435         return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
0436     }
0437
0438     const deriv_type& get_current_deriv( void ) const
0439     {
0440         return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
0441     }
0442
0443     deriv_type& get_old_deriv( void )
0444     {
0445         return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
0446     }
0447
0448     const deriv_type& get_old_deriv( void ) const
0449     {
0450         return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
0451     }
0452
0453
0454     void toggle_current_state( void )
0455     {
0456         m_current_state_x1 = ! m_current_state_x1;
0457     }
0458
0459
0460     controlled_stepper_type m_stepper;
0461     resizer_type m_resizer;
0462     bool m_current_state_x1;
0463     wrapped_state_type m_x1 , m_x2;
0464     wrapped_deriv_type m_dxdt1 , m_dxdt2;
0465     time_type m_t , m_t_old , m_dt;
0466     bool m_is_deriv_initialized;
0467
0468 };
0469
0470 } // namespace odeint
0471 } // namespace numeric
0472 } // namespace boost
0473
0474
0475
0476 #endif // BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED