Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:42:56 

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     BOOST_STATIC_ASSERT_MSG(
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(
0176             in, detail::bind( &stepper_type::template resize_impl< StateIn >,
0177                               detail::ref( *this ), detail::_1 ) );
0178         size_t k = 0;
0179         m_midpoint.set_steps( m_interval_sequence[k] );
0180         m_midpoint.do_step( system, in, dxdt, t, out, dt );
0181         for ( k = 1; k <= m_k_max; ++k )
0182         {
0183             m_midpoint.set_steps( m_interval_sequence[k] );
0184             m_midpoint.do_step( system, in, dxdt, t, m_table[k - 1].m_v, dt );
0185             extrapolate( k, m_table, m_coeff, out );
0186         }
0187     }
0188 
0189     template < class System, class StateInOut, class DerivIn >
0190     void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt,
0191                           time_type t, time_type dt )
0192     {
0193         // special care for inout
0194         m_xout_resizer.adjust_size(
0195             inout,
0196             detail::bind( &stepper_type::template resize_m_xout< StateInOut >,
0197                           detail::ref( *this ), detail::_1 ) );
0198         do_step_impl( system, inout, dxdt, t, m_xout.m_v, dt );
0199         boost::numeric::odeint::copy( m_xout.m_v, inout );
0200     }
0201 
0202     template < class System, class StateInOut, class DerivIn >
0203     void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt,
0204                             time_type t, time_type dt )
0205     {
0206         do_step_impl_io( system , x , dxdt , t , dt );
0207     }
0208 
0209     template < class System, class StateIn, class DerivIn, class StateOut >
0210     void do_step_dxdt_impl( System system, const StateIn &in,
0211                             const DerivIn &dxdt, time_type t, StateOut &out,
0212                             time_type dt )
0213     {
0214         do_step_impl( system , in , dxdt , t , out , dt );
0215     }
0216 
0217 
0218     template < class StateIn > void adjust_size( const StateIn &x )
0219     {
0220         resize_impl( x );
0221         m_midpoint.adjust_size( x );
0222     }
0223 
0224   private:
0225     template < class StateIn > bool resize_impl( const StateIn &x )
0226     {
0227         bool resized( false );
0228         for ( size_t i = 0; i < m_k_max; ++i )
0229             resized |= adjust_size_by_resizeability(
0230                 m_table[i], x, typename is_resizeable< state_type >::type() );
0231         return resized;
0232     }
0233 
0234     template < class StateIn > bool resize_m_xout( const StateIn &x )
0235     {
0236         return adjust_size_by_resizeability(
0237             m_xout, x, typename is_resizeable< state_type >::type() );
0238     }
0239 
0240     template < class StateInOut >
0241     void extrapolate( size_t k, state_table_type &table,
0242                       const value_matrix &coeff, StateInOut &xest )
0243     /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
0244        uses the obtained intermediate results to extrapolate to dt->0
0245     */
0246     {
0247         static const value_type val1 = static_cast< value_type >( 1.0 );
0248 
0249         for ( int j = k - 1; j > 0; --j )
0250         {
0251             stepper_base_type::m_algebra.for_each3(
0252                 table[j - 1].m_v, table[j].m_v, table[j - 1].m_v,
0253                 typename operations_type::template scale_sum2<
0254                     value_type, value_type >( val1 + coeff[k][j],
0255                                               -coeff[k][j] ) );
0256         }
0257         stepper_base_type::m_algebra.for_each3(
0258             xest, table[0].m_v, xest,
0259             typename operations_type::template scale_sum2<
0260                 value_type, value_type >( val1 + coeff[k][0], -coeff[k][0] ) );
0261     }
0262 
0263   private:
0264     midpoint_stepper_type m_midpoint;
0265 
0266     resizer_type m_resizer;
0267     resizer_type m_xout_resizer;
0268 
0269     int_vector m_interval_sequence; // stores the successive interval counts
0270     value_matrix m_coeff;
0271 
0272     wrapped_state_type m_xout;
0273     state_table_type m_table; // sequence of states for extrapolation
0274 };
0275 
0276 /******** DOXYGEN *******/
0277 
0278 /**
0279  * \class extrapolation_stepper
0280  * \brief Extrapolation stepper with configurable order, and error estimation.
0281  *
0282  * The extrapolation stepper is a stepper with error estimation and configurable
0283  * order. The order is given as template parameter and needs to be an _odd_
0284  * number. The stepper is based on several executions of the modified midpoint
0285  * method and a Richardson extrapolation. This is essentially the same technique
0286  * as for bulirsch_stoer, but without the variable order.
0287  *
0288  * \note The Order parameter has to be an even number greater 2.
0289  */
0290 }
0291 }
0292 }
0293 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!