Chris@102: /*
Chris@102:   [auto_generated]
Chris@102:   boost/numeric/odeint/stepper/extrapolation_stepper.hpp
Chris@102: 
Chris@102:   [begin_description]
Chris@102:   extrapolation stepper
Chris@102:   [end_description]
Chris@102: 
Chris@102:   Copyright 2009-2015 Mario Mulansky
Chris@102: 
Chris@102:   Distributed under the Boost Software License, Version 1.0.
Chris@102:   (See accompanying file LICENSE_1_0.txt or
Chris@102:   copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@102: */
Chris@102: 
Chris@102: #ifndef BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED
Chris@102: #define BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED
Chris@102: 
Chris@102: #include <iostream>
Chris@102: 
Chris@102: #include <algorithm>
Chris@102: 
Chris@102: #include <boost/config.hpp> // for min/max guidelines
Chris@102: #include <boost/static_assert.hpp>
Chris@102: 
Chris@102: #include <boost/numeric/odeint/util/bind.hpp>
Chris@102: #include <boost/numeric/odeint/util/unwrap_reference.hpp>
Chris@102: 
Chris@102: #include <boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp>
Chris@102: #include <boost/numeric/odeint/stepper/modified_midpoint.hpp>
Chris@102: #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
Chris@102: #include <boost/numeric/odeint/algebra/range_algebra.hpp>
Chris@102: #include <boost/numeric/odeint/algebra/default_operations.hpp>
Chris@102: #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
Chris@102: #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
Chris@102: 
Chris@102: #include <boost/numeric/odeint/util/state_wrapper.hpp>
Chris@102: #include <boost/numeric/odeint/util/is_resizeable.hpp>
Chris@102: #include <boost/numeric/odeint/util/resizer.hpp>
Chris@102: #include <boost/numeric/odeint/util/unit_helper.hpp>
Chris@102: #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
Chris@102: 
Chris@102: namespace boost
Chris@102: {
Chris@102: namespace numeric
Chris@102: {
Chris@102: namespace odeint
Chris@102: {
Chris@102: 
Chris@102: template < unsigned short Order, class State, class Value = double,
Chris@102:            class Deriv = State, class Time = Value,
Chris@102:            class Algebra = typename algebra_dispatcher< State >::algebra_type,
Chris@102:            class Operations =
Chris@102:                typename operations_dispatcher< State >::operations_type,
Chris@102:            class Resizer = initially_resizer >
Chris@102: #ifndef DOXYGEN_SKIP
Chris@102: class extrapolation_stepper
Chris@102:     : public explicit_error_stepper_base<
Chris@102:           extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra,
Chris@102:                                  Operations, Resizer >,
Chris@102:           Order, Order, Order - 2, State, Value, Deriv, Time, Algebra,
Chris@102:           Operations, Resizer >
Chris@102: #else
Chris@102: class extrapolation_stepper : public explicit_error_stepper_base
Chris@102: #endif
Chris@102: {
Chris@102: 
Chris@102:   private:
Chris@102:     // check for Order being odd
Chris@102:     BOOST_STATIC_ASSERT_MSG(
Chris@102:         ( ( Order % 2 ) == 0 ) && ( Order > 2 ),
Chris@102:         "extrapolation_stepper requires even Order larger than 2" );
Chris@102: 
Chris@102:   public:
Chris@102: #ifndef DOXYGEN_SKIP
Chris@102:     typedef explicit_error_stepper_base<
Chris@102:         extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra,
Chris@102:                                Operations, Resizer >,
Chris@102:         Order, Order, Order - 2, State, Value, Deriv, Time, Algebra, Operations,
Chris@102:         Resizer > stepper_base_type;
Chris@102: #else
Chris@102:     typedef explicit_error_stepper_base< extrapolation_stepper< ... >, ... >
Chris@102:     stepper_base_type;
Chris@102: #endif
Chris@102: 
Chris@102:     typedef typename stepper_base_type::state_type state_type;
Chris@102:     typedef typename stepper_base_type::value_type value_type;
Chris@102:     typedef typename stepper_base_type::deriv_type deriv_type;
Chris@102:     typedef typename stepper_base_type::time_type time_type;
Chris@102:     typedef typename stepper_base_type::algebra_type algebra_type;
Chris@102:     typedef typename stepper_base_type::operations_type operations_type;
Chris@102:     typedef typename stepper_base_type::resizer_type resizer_type;
Chris@102: 
Chris@102: #ifndef DOXYGEN_SKIP
Chris@102:     typedef typename stepper_base_type::stepper_type stepper_type;
Chris@102:     typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
Chris@102:     typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
Chris@102: 
Chris@102:     typedef std::vector< value_type > value_vector;
Chris@102:     typedef std::vector< value_vector > value_matrix;
Chris@102:     typedef std::vector< size_t > int_vector;
Chris@102:     typedef std::vector< wrapped_state_type > state_table_type;
Chris@102:     typedef modified_midpoint< state_type, value_type, deriv_type, time_type,
Chris@102:                                algebra_type, operations_type,
Chris@102:                                resizer_type > midpoint_stepper_type;
Chris@102: 
Chris@102: #endif // DOXYGEN_SKIP
Chris@102: 
Chris@102:     typedef unsigned short order_type;
Chris@102:     static const order_type order_value = stepper_base_type::order_value;
Chris@102:     static const order_type stepper_order_value =
Chris@102:         stepper_base_type::stepper_order_value;
Chris@102:     static const order_type error_order_value =
Chris@102:         stepper_base_type::error_order_value;
Chris@102: 
Chris@102:     const static size_t m_k_max = ( order_value - 2 ) / 2;
Chris@102: 
Chris@102:     extrapolation_stepper( const algebra_type &algebra = algebra_type() )
Chris@102:         : stepper_base_type( algebra ), m_interval_sequence( m_k_max + 1 ),
Chris@102:           m_coeff( m_k_max + 1 ), m_table( m_k_max )
Chris@102:     {
Chris@102:         for ( unsigned short i = 0; i < m_k_max + 1; i++ )
Chris@102:         {
Chris@102:             m_interval_sequence[i] = 2 * ( i + 1 );
Chris@102:             m_coeff[i].resize( i );
Chris@102:             for ( size_t k = 0; k < i; ++k )
Chris@102:             {
Chris@102:                 const value_type r =
Chris@102:                     static_cast< value_type >( m_interval_sequence[i] ) /
Chris@102:                     static_cast< value_type >( m_interval_sequence[k] );
Chris@102:                 m_coeff[i][k] =
Chris@102:                     static_cast< value_type >( 1 ) /
Chris@102:                     ( r * r - static_cast< value_type >(
Chris@102:                                   1 ) ); // coefficients for extrapolation
Chris@102:             }
Chris@102:         }
Chris@102:     }
Chris@102: 
Chris@102:     template < class System, class StateIn, class DerivIn, class StateOut,
Chris@102:                class Err >
Chris@102:     void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt,
Chris@102:                        time_type t, StateOut &out, time_type dt, Err &xerr )
Chris@102:     {
Chris@102:         // std::cout << "dt: " << dt << std::endl;
Chris@102:         // normal step
Chris@102:         do_step_impl( system, in, dxdt, t, out, dt );
Chris@102: 
Chris@102:         static const value_type val1( 1.0 );
Chris@102:         // additionally, perform the error calculation
Chris@102:         stepper_base_type::m_algebra.for_each3(
Chris@102:             xerr, out, m_table[0].m_v,
Chris@102:             typename operations_type::template scale_sum2<
Chris@102:                 value_type, value_type >( val1, -val1 ) );
Chris@102:     }
Chris@102: 
Chris@102:     template < class System, class StateInOut, class DerivIn, class Err >
Chris@102:     void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt,
Chris@102:                           time_type t, time_type dt, Err &xerr )
Chris@102:     {
Chris@102:         // normal step
Chris@102:         do_step_impl_io( system, inout, dxdt, t, dt );
Chris@102: 
Chris@102:         static const value_type val1( 1.0 );
Chris@102:         // additionally, perform the error calculation
Chris@102:         stepper_base_type::m_algebra.for_each3(
Chris@102:             xerr, inout, m_table[0].m_v,
Chris@102:             typename operations_type::template scale_sum2<
Chris@102:                 value_type, value_type >( val1, -val1 ) );
Chris@102:     }
Chris@102: 
Chris@102:     template < class System, class StateIn, class DerivIn, class StateOut >
Chris@102:     void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt,
Chris@102:                        time_type t, StateOut &out, time_type dt )
Chris@102:     {
Chris@102:         m_resizer.adjust_size(
Chris@102:             in, detail::bind( &stepper_type::template resize_impl< StateIn >,
Chris@102:                               detail::ref( *this ), detail::_1 ) );
Chris@102:         size_t k = 0;
Chris@102:         m_midpoint.set_steps( m_interval_sequence[k] );
Chris@102:         m_midpoint.do_step( system, in, dxdt, t, out, dt );
Chris@102:         for ( k = 1; k <= m_k_max; ++k )
Chris@102:         {
Chris@102:             m_midpoint.set_steps( m_interval_sequence[k] );
Chris@102:             m_midpoint.do_step( system, in, dxdt, t, m_table[k - 1].m_v, dt );
Chris@102:             extrapolate( k, m_table, m_coeff, out );
Chris@102:         }
Chris@102:     }
Chris@102: 
Chris@102:     template < class System, class StateInOut, class DerivIn >
Chris@102:     void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt,
Chris@102:                           time_type t, time_type dt )
Chris@102:     {
Chris@102:         // special care for inout
Chris@102:         m_xout_resizer.adjust_size(
Chris@102:             inout,
Chris@102:             detail::bind( &stepper_type::template resize_m_xout< StateInOut >,
Chris@102:                           detail::ref( *this ), detail::_1 ) );
Chris@102:         do_step_impl( system, inout, dxdt, t, m_xout.m_v, dt );
Chris@102:         boost::numeric::odeint::copy( m_xout.m_v, inout );
Chris@102:     }
Chris@102: 
Chris@102:     template < class System, class StateInOut, class DerivIn >
Chris@102:     void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt,
Chris@102:                             time_type t, time_type dt )
Chris@102:     {
Chris@102:         do_step_impl_io( system , x , dxdt , t , dt );
Chris@102:     }
Chris@102: 
Chris@102:     template < class System, class StateIn, class DerivIn, class StateOut >
Chris@102:     void do_step_dxdt_impl( System system, const StateIn &in,
Chris@102:                             const DerivIn &dxdt, time_type t, StateOut &out,
Chris@102:                             time_type dt )
Chris@102:     {
Chris@102:         do_step_impl( system , in , dxdt , t , out , dt );
Chris@102:     }
Chris@102: 
Chris@102: 
Chris@102:     template < class StateIn > void adjust_size( const StateIn &x )
Chris@102:     {
Chris@102:         resize_impl( x );
Chris@102:         m_midpoint.adjust_size( x );
Chris@102:     }
Chris@102: 
Chris@102:   private:
Chris@102:     template < class StateIn > bool resize_impl( const StateIn &x )
Chris@102:     {
Chris@102:         bool resized( false );
Chris@102:         for ( size_t i = 0; i < m_k_max; ++i )
Chris@102:             resized |= adjust_size_by_resizeability(
Chris@102:                 m_table[i], x, typename is_resizeable< state_type >::type() );
Chris@102:         return resized;
Chris@102:     }
Chris@102: 
Chris@102:     template < class StateIn > bool resize_m_xout( const StateIn &x )
Chris@102:     {
Chris@102:         return adjust_size_by_resizeability(
Chris@102:             m_xout, x, typename is_resizeable< state_type >::type() );
Chris@102:     }
Chris@102: 
Chris@102:     template < class StateInOut >
Chris@102:     void extrapolate( size_t k, state_table_type &table,
Chris@102:                       const value_matrix &coeff, StateInOut &xest )
Chris@102:     /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
Chris@102:        uses the obtained intermediate results to extrapolate to dt->0
Chris@102:     */
Chris@102:     {
Chris@102:         static const value_type val1 = static_cast< value_type >( 1.0 );
Chris@102: 
Chris@102:         for ( int j = k - 1; j > 0; --j )
Chris@102:         {
Chris@102:             stepper_base_type::m_algebra.for_each3(
Chris@102:                 table[j - 1].m_v, table[j].m_v, table[j - 1].m_v,
Chris@102:                 typename operations_type::template scale_sum2<
Chris@102:                     value_type, value_type >( val1 + coeff[k][j],
Chris@102:                                               -coeff[k][j] ) );
Chris@102:         }
Chris@102:         stepper_base_type::m_algebra.for_each3(
Chris@102:             xest, table[0].m_v, xest,
Chris@102:             typename operations_type::template scale_sum2<
Chris@102:                 value_type, value_type >( val1 + coeff[k][0], -coeff[k][0] ) );
Chris@102:     }
Chris@102: 
Chris@102:   private:
Chris@102:     midpoint_stepper_type m_midpoint;
Chris@102: 
Chris@102:     resizer_type m_resizer;
Chris@102:     resizer_type m_xout_resizer;
Chris@102: 
Chris@102:     int_vector m_interval_sequence; // stores the successive interval counts
Chris@102:     value_matrix m_coeff;
Chris@102: 
Chris@102:     wrapped_state_type m_xout;
Chris@102:     state_table_type m_table; // sequence of states for extrapolation
Chris@102: };
Chris@102: 
Chris@102: /******** DOXYGEN *******/
Chris@102: 
Chris@102: /**
Chris@102:  * \class extrapolation_stepper
Chris@102:  * \brief Extrapolation stepper with configurable order, and error estimation.
Chris@102:  *
Chris@102:  * The extrapolation stepper is a stepper with error estimation and configurable
Chris@102:  * order. The order is given as template parameter and needs to be an _odd_
Chris@102:  * number. The stepper is based on several executions of the modified midpoint
Chris@102:  * method and a Richardson extrapolation. This is essentially the same technique
Chris@102:  * as for bulirsch_stoer, but without the variable order.
Chris@102:  *
Chris@102:  * \note The Order parameter has to be an even number greater 2.
Chris@102:  */
Chris@102: }
Chris@102: }
Chris@102: }
Chris@102: #endif