odeint/stepper/extrapolation_stepper.hpp

0001 /*
0002   [auto_generated]
0003   boost/numeric/odeint/stepper/extrapolation_stepper.hpp
0004
0005   [begin_description]
0006   extrapolation stepper
0007   [end_description]
0008
0009   Copyright 2009-2015 Mario Mulansky
0010
0011   Distributed under the Boost Software License, Version 1.0.
0012   (See accompanying file LICENSE_1_0.txt or
0013   copy at http://www.boost.org/LICENSE_1_0.txt)
0014 */
0015
0016 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED
0017 #define BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED
0018
0019 #include <iostream>
0020
0021 #include <algorithm>
0022
0023 #include <boost/config.hpp> // for min/max guidelines
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/stepper/base/explicit_error_stepper_base.hpp>
0030 #include <boost/numeric/odeint/stepper/modified_midpoint.hpp>
0031 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
0032 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
0033 #include <boost/numeric/odeint/algebra/default_operations.hpp>
0034 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
0035 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
0036
0037 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0038 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0039 #include <boost/numeric/odeint/util/resizer.hpp>
0040 #include <boost/numeric/odeint/util/unit_helper.hpp>
0041 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
0042
0043 namespace boost
0044 {
0045 namespace numeric
0046 {
0047 namespace odeint
0048 {
0049
0050 template < unsigned short Order, class State, class Value = double,
0051            class Deriv = State, class Time = Value,
0052            class Algebra = typename algebra_dispatcher< State >::algebra_type,
0053            class Operations =
0054                typename operations_dispatcher< State >::operations_type,
0055            class Resizer = initially_resizer >
0056 #ifndef DOXYGEN_SKIP
0057 class extrapolation_stepper
0058     : public explicit_error_stepper_base<
0059           extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra,
0060                                  Operations, Resizer >,
0061           Order, Order, Order - 2, State, Value, Deriv, Time, Algebra,
0062           Operations, Resizer >
0063 #else
0064 class extrapolation_stepper : public explicit_error_stepper_base
0065 #endif
0066 {
0067
0068   private:
0069     // check for Order being odd
0070     static_assert(
0071         ( ( Order % 2 ) == 0 ) && ( Order > 2 ),
0072         "extrapolation_stepper requires even Order larger than 2" );
0073
0074   public:
0075 #ifndef DOXYGEN_SKIP
0076     typedef explicit_error_stepper_base<
0077         extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra,
0078                                Operations, Resizer >,
0079         Order, Order, Order - 2, State, Value, Deriv, Time, Algebra, Operations,
0080         Resizer > stepper_base_type;
0081 #else
0082     typedef explicit_error_stepper_base< extrapolation_stepper< ... >, ... >
0083     stepper_base_type;
0084 #endif
0085
0086     typedef typename stepper_base_type::state_type state_type;
0087     typedef typename stepper_base_type::value_type value_type;
0088     typedef typename stepper_base_type::deriv_type deriv_type;
0089     typedef typename stepper_base_type::time_type time_type;
0090     typedef typename stepper_base_type::algebra_type algebra_type;
0091     typedef typename stepper_base_type::operations_type operations_type;
0092     typedef typename stepper_base_type::resizer_type resizer_type;
0093
0094 #ifndef DOXYGEN_SKIP
0095     typedef typename stepper_base_type::stepper_type stepper_type;
0096     typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
0097     typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
0098
0099     typedef std::vector< value_type > value_vector;
0100     typedef std::vector< value_vector > value_matrix;
0101     typedef std::vector< size_t > int_vector;
0102     typedef std::vector< wrapped_state_type > state_table_type;
0103     typedef modified_midpoint< state_type, value_type, deriv_type, time_type,
0104                                algebra_type, operations_type,
0105                                resizer_type > midpoint_stepper_type;
0106
0107 #endif // DOXYGEN_SKIP
0108
0109     typedef unsigned short order_type;
0110     static const order_type order_value = stepper_base_type::order_value;
0111     static const order_type stepper_order_value =
0112         stepper_base_type::stepper_order_value;
0113     static const order_type error_order_value =
0114         stepper_base_type::error_order_value;
0115
0116     const static size_t m_k_max = ( order_value - 2 ) / 2;
0117
0118     extrapolation_stepper( const algebra_type &algebra = algebra_type() )
0119         : stepper_base_type( algebra ), m_interval_sequence( m_k_max + 1 ),
0120           m_coeff( m_k_max + 1 ), m_table( m_k_max )
0121     {
0122         for ( unsigned short i = 0; i < m_k_max + 1; i++ )
0123         {
0124             m_interval_sequence[i] = 2 * ( i + 1 );
0125             m_coeff[i].resize( i );
0126             for ( size_t k = 0; k < i; ++k )
0127             {
0128                 const value_type r =
0129                     static_cast< value_type >( m_interval_sequence[i] ) /
0130                     static_cast< value_type >( m_interval_sequence[k] );
0131                 m_coeff[i][k] =
0132                     static_cast< value_type >( 1 ) /
0133                     ( r * r - static_cast< value_type >(
0134                                   1 ) ); // coefficients for extrapolation
0135             }
0136         }
0137     }
0138
0139     template < class System, class StateIn, class DerivIn, class StateOut,
0140                class Err >
0141     void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt,
0142                        time_type t, StateOut &out, time_type dt, Err &xerr )
0143     {
0144         // std::cout << "dt: " << dt << std::endl;
0145         // normal step
0146         do_step_impl( system, in, dxdt, t, out, dt );
0147
0148         static const value_type val1( 1.0 );
0149         // additionally, perform the error calculation
0150         stepper_base_type::m_algebra.for_each3(
0151             xerr, out, m_table[0].m_v,
0152             typename operations_type::template scale_sum2<
0153                 value_type, value_type >( val1, -val1 ) );
0154     }
0155
0156     template < class System, class StateInOut, class DerivIn, class Err >
0157     void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt,
0158                           time_type t, time_type dt, Err &xerr )
0159     {
0160         // normal step
0161         do_step_impl_io( system, inout, dxdt, t, dt );
0162
0163         static const value_type val1( 1.0 );
0164         // additionally, perform the error calculation
0165         stepper_base_type::m_algebra.for_each3(
0166             xerr, inout, m_table[0].m_v,
0167             typename operations_type::template scale_sum2<
0168                 value_type, value_type >( val1, -val1 ) );
0169     }
0170
0171     template < class System, class StateIn, class DerivIn, class StateOut >
0172     void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt,
0173                        time_type t, StateOut &out, time_type dt )
0174     {
0175         m_resizer.adjust_size(in, [this](auto&& arg) { return this->resize_impl<StateIn>(std::forward<decltype(arg)>(arg)); });
0176         size_t k = 0;
0177         m_midpoint.set_steps( m_interval_sequence[k] );
0178         m_midpoint.do_step( system, in, dxdt, t, out, dt );
0179         for ( k = 1; k <= m_k_max; ++k )
0180         {
0181             m_midpoint.set_steps( m_interval_sequence[k] );
0182             m_midpoint.do_step( system, in, dxdt, t, m_table[k - 1].m_v, dt );
0183             extrapolate( k, m_table, m_coeff, out );
0184         }
0185     }
0186
0187     template < class System, class StateInOut, class DerivIn >
0188     void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt,
0189                           time_type t, time_type dt )
0190     {
0191         // special care for inout
0192         m_xout_resizer.adjust_size(inout, [this](auto&& arg) { return this->resize_m_xout<StateInOut>(std::forward<decltype(arg)>(arg)); });
0193         do_step_impl( system, inout, dxdt, t, m_xout.m_v, dt );
0194         boost::numeric::odeint::copy( m_xout.m_v, inout );
0195     }
0196
0197     template < class System, class StateInOut, class DerivIn >
0198     void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt,
0199                             time_type t, time_type dt )
0200     {
0201         do_step_impl_io( system , x , dxdt , t , dt );
0202     }
0203
0204     template < class System, class StateIn, class DerivIn, class StateOut >
0205     void do_step_dxdt_impl( System system, const StateIn &in,
0206                             const DerivIn &dxdt, time_type t, StateOut &out,
0207                             time_type dt )
0208     {
0209         do_step_impl( system , in , dxdt , t , out , dt );
0210     }
0211
0212
0213     template < class StateIn > void adjust_size( const StateIn &x )
0214     {
0215         resize_impl( x );
0216         m_midpoint.adjust_size( x );
0217     }
0218
0219   private:
0220     template < class StateIn > bool resize_impl( const StateIn &x )
0221     {
0222         bool resized( false );
0223         for ( size_t i = 0; i < m_k_max; ++i )
0224             resized |= adjust_size_by_resizeability(
0225                 m_table[i], x, typename is_resizeable< state_type >::type() );
0226         return resized;
0227     }
0228
0229     template < class StateIn > bool resize_m_xout( const StateIn &x )
0230     {
0231         return adjust_size_by_resizeability(
0232             m_xout, x, typename is_resizeable< state_type >::type() );
0233     }
0234
0235     template < class StateInOut >
0236     void extrapolate( size_t k, state_table_type &table,
0237                       const value_matrix &coeff, StateInOut &xest )
0238     /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
0239        uses the obtained intermediate results to extrapolate to dt->0
0240     */
0241     {
0242         static const value_type val1 = static_cast< value_type >( 1.0 );
0243
0244         for ( int j = k - 1; j > 0; --j )
0245         {
0246             stepper_base_type::m_algebra.for_each3(
0247                 table[j - 1].m_v, table[j].m_v, table[j - 1].m_v,
0248                 typename operations_type::template scale_sum2<
0249                     value_type, value_type >( val1 + coeff[k][j],
0250                                               -coeff[k][j] ) );
0251         }
0252         stepper_base_type::m_algebra.for_each3(
0253             xest, table[0].m_v, xest,
0254             typename operations_type::template scale_sum2<
0255                 value_type, value_type >( val1 + coeff[k][0], -coeff[k][0] ) );
0256     }
0257
0258   private:
0259     midpoint_stepper_type m_midpoint;
0260
0261     resizer_type m_resizer;
0262     resizer_type m_xout_resizer;
0263
0264     int_vector m_interval_sequence; // stores the successive interval counts
0265     value_matrix m_coeff;
0266
0267     wrapped_state_type m_xout;
0268     state_table_type m_table; // sequence of states for extrapolation
0269 };
0270
0271 /******** DOXYGEN *******/
0272
0273 /**
0274  * \class extrapolation_stepper
0275  * \brief Extrapolation stepper with configurable order, and error estimation.
0276  *
0277  * The extrapolation stepper is a stepper with error estimation and configurable
0278  * order. The order is given as template parameter and needs to be an _odd_
0279  * number. The stepper is based on several executions of the modified midpoint
0280  * method and a Richardson extrapolation. This is essentially the same technique
0281  * as for bulirsch_stoer, but without the variable order.
0282  *
0283  * \note The Order parameter has to be an even number greater 2.
0284  */
0285 }
0286 }
0287 }
0288 #endif