odeint/stepper/bulirsch_stoer_dense_out.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
0004
0005  [begin_description]
0006  Implementaiton of the Burlish-Stoer method with dense output
0007  [end_description]
0008
0009  Copyright 2011-2015 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_BULIRSCH_STOER_DENSE_OUT_HPP_INCLUDED
0020 #define BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_DENSE_OUT_HPP_INCLUDED
0021
0022
0023 #include <iostream>
0024
0025 #include <algorithm>
0026
0027 #include <boost/config.hpp> // for min/max guidelines
0028
0029 #include <boost/numeric/odeint/util/bind.hpp>
0030
0031 #include <boost/math/special_functions/binomial.hpp>
0032
0033 #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
0034 #include <boost/numeric/odeint/stepper/modified_midpoint.hpp>
0035 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
0036 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
0037 #include <boost/numeric/odeint/algebra/default_operations.hpp>
0038 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
0039 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
0040
0041 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0042 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0043 #include <boost/numeric/odeint/util/resizer.hpp>
0044 #include <boost/numeric/odeint/util/unit_helper.hpp>
0045
0046 #include <boost/numeric/odeint/integrate/max_step_checker.hpp>
0047
0048 #include <boost/type_traits.hpp>
0049
0050
0051 namespace boost {
0052 namespace numeric {
0053 namespace odeint {
0054
0055 template<
0056     class State ,
0057     class Value = double ,
0058     class Deriv = State ,
0059     class Time = Value ,
0060     class Algebra = typename algebra_dispatcher< State >::algebra_type ,
0061     class Operations = typename operations_dispatcher< State >::operations_type ,
0062     class Resizer = initially_resizer
0063     >
0064 class bulirsch_stoer_dense_out {
0065
0066
0067 public:
0068
0069     typedef State state_type;
0070     typedef Value value_type;
0071     typedef Deriv deriv_type;
0072     typedef Time time_type;
0073     typedef Algebra algebra_type;
0074     typedef Operations operations_type;
0075     typedef Resizer resizer_type;
0076     typedef dense_output_stepper_tag stepper_category;
0077 #ifndef DOXYGEN_SKIP
0078     typedef state_wrapper< state_type > wrapped_state_type;
0079     typedef state_wrapper< deriv_type > wrapped_deriv_type;
0080
0081     typedef bulirsch_stoer_dense_out< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type;
0082
0083     typedef typename inverse_time< time_type >::type inv_time_type;
0084
0085     typedef std::vector< value_type > value_vector;
0086     typedef std::vector< time_type > time_vector;
0087     typedef std::vector< inv_time_type > inv_time_vector;  //should be 1/time_type for boost.units
0088     typedef std::vector< value_vector > value_matrix;
0089     typedef std::vector< size_t > int_vector;
0090     typedef std::vector< wrapped_state_type > state_vector_type;
0091     typedef std::vector< wrapped_deriv_type > deriv_vector_type;
0092     typedef std::vector< deriv_vector_type > deriv_table_type;
0093 #endif //DOXYGEN_SKIP
0094
0095     const static size_t m_k_max = 8;
0096
0097
0098
0099     bulirsch_stoer_dense_out(
0100         value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
0101         value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 ,
0102         time_type max_dt = static_cast<time_type>(0) ,
0103         bool control_interpolation = false )
0104         : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) ,
0105           m_max_dt(max_dt) ,
0106           m_control_interpolation( control_interpolation) ,
0107           m_last_step_rejected( false ) , m_first( true ) ,
0108           m_current_state_x1( true ) ,
0109           m_error( m_k_max ) ,
0110           m_interval_sequence( m_k_max+1 ) ,
0111           m_coeff( m_k_max+1 ) ,
0112           m_cost( m_k_max+1 ) ,
0113           m_facmin_table( m_k_max+1 ) ,
0114           m_table( m_k_max ) ,
0115           m_mp_states( m_k_max+1 ) ,
0116           m_derivs( m_k_max+1 ) ,
0117           m_diffs( 2*m_k_max+2 ) ,
0118           STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 )
0119     {
0120         BOOST_USING_STD_MIN();
0121         BOOST_USING_STD_MAX();
0122
0123         for( unsigned short i = 0; i < m_k_max+1; i++ )
0124         {
0125             /* only this specific sequence allows for dense output */
0126             m_interval_sequence[i] = 2 + 4*i;  // 2 6 10 14 ...
0127             m_derivs[i].resize( m_interval_sequence[i] );
0128             if( i == 0 )
0129             {
0130                 m_cost[i] = m_interval_sequence[i];
0131             } else
0132             {
0133                 m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
0134             }
0135             m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) );
0136             m_coeff[i].resize(i);
0137             for( size_t k = 0 ; k < i ; ++k  )
0138             {
0139                 const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
0140                 m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
0141             }
0142             // crude estimate of optimal order
0143
0144             m_current_k_opt = 4;
0145             /* no calculation because log10 might not exist for value_type!
0146             const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >( 1.0E-12 ) ) ) * 0.6 + 0.5 );
0147             m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 1 , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>( m_k_max-1 ) , static_cast<int>( logfact ) ));
0148             */
0149         }
0150         int num = 1;
0151         for( int i = 2*(m_k_max)+1 ; i >=0  ; i-- )
0152         {
0153             m_diffs[i].resize( num );
0154             num += (i+1)%2;
0155         }
0156     }
0157
0158     template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
0159     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , DerivOut &dxdt_new , time_type &dt )
0160     {
0161         if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
0162         {
0163             // given step size is bigger then max_dt
0164             // set limit and return fail
0165             dt = m_max_dt;
0166             return fail;
0167         }
0168
0169         BOOST_USING_STD_MIN();
0170         BOOST_USING_STD_MAX();
0171         using std::pow;
0172
0173         static const value_type val1( 1.0 );
0174
0175         bool reject( true );
0176
0177         time_vector h_opt( m_k_max+1 );
0178         inv_time_vector work( m_k_max+1 );
0179
0180         m_k_final = 0;
0181         time_type new_h = dt;
0182
0183         //std::cout << "t=" << t <<", dt=" << dt << ", k_opt=" << m_current_k_opt << ", first: " << m_first << std::endl;
0184
0185         for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ )
0186         {
0187             m_midpoint.set_steps( m_interval_sequence[k] );
0188             if( k == 0 )
0189             {
0190                 m_midpoint.do_step( system , in , dxdt , t , out , dt , m_mp_states[k].m_v , m_derivs[k]);
0191             }
0192             else
0193             {
0194                 m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt , m_mp_states[k].m_v , m_derivs[k] );
0195                 extrapolate( k , m_table , m_coeff , out );
0196                 // get error estimate
0197                 m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v ,
0198                                      typename operations_type::template scale_sum2< value_type , value_type >( val1 , -val1 ) );
0199                 const value_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt );
0200                 h_opt[k] = calc_h_opt( dt , error , k );
0201                 work[k] = static_cast<value_type>( m_cost[k] ) / h_opt[k];
0202
0203                 m_k_final = k;
0204
0205                 if( (k == m_current_k_opt-1) || m_first )
0206                 { // convergence before k_opt ?
0207                     if( error < 1.0 )
0208                     {
0209                         //convergence
0210                         reject = false;
0211                         if( (work[k] < KFAC2*work[k-1]) || (m_current_k_opt <= 2) )
0212                         {
0213                             // leave order as is (except we were in first round)
0214                             m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(k)+1 ) );
0215                             new_h = h_opt[k] * static_cast<value_type>( m_cost[k+1] ) / static_cast<value_type>( m_cost[k] );
0216                         } else {
0217                             m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(k) ) );
0218                             new_h = h_opt[k];
0219                         }
0220                         break;
0221                     }
0222                     else if( should_reject( error , k ) && !m_first )
0223                     {
0224                         reject = true;
0225                         new_h = h_opt[k];
0226                         break;
0227                     }
0228                 }
0229                 if( k == m_current_k_opt )
0230                 { // convergence at k_opt ?
0231                     if( error < 1.0 )
0232                     {
0233                         //convergence
0234                         reject = false;
0235                         if( (work[k-1] < KFAC2*work[k]) )
0236                         {
0237                             m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
0238                             new_h = h_opt[m_current_k_opt];
0239                         }
0240                         else if( (work[k] < KFAC2*work[k-1]) && !m_last_step_rejected )
0241                         {
0242                             m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , static_cast<int>(m_current_k_opt)+1 );
0243                             new_h = h_opt[k]*static_cast<value_type>( m_cost[m_current_k_opt] ) / static_cast<value_type>( m_cost[k] );
0244                         } else
0245                             new_h = h_opt[m_current_k_opt];
0246                         break;
0247                     }
0248                     else if( should_reject( error , k ) )
0249                     {
0250                         reject = true;
0251                         new_h = h_opt[m_current_k_opt];
0252                         break;
0253                     }
0254                 }
0255                 if( k == m_current_k_opt+1 )
0256                 { // convergence at k_opt+1 ?
0257                     if( error < 1.0 )
0258                     {   //convergence
0259                         reject = false;
0260                         if( work[k-2] < KFAC2*work[k-1] )
0261                             m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
0262                         if( (work[k] < KFAC2*work[m_current_k_opt]) && !m_last_step_rejected )
0263                             m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , static_cast<int>(k) );
0264                         new_h = h_opt[m_current_k_opt];
0265                     } else
0266                     {
0267                         reject = true;
0268                         new_h = h_opt[m_current_k_opt];
0269                     }
0270                     break;
0271                 }
0272             }
0273         }
0274
0275         if( !reject )
0276         {
0277
0278             //calculate dxdt for next step and dense output
0279             typename odeint::unwrap_reference< System >::type &sys = system;
0280             sys( out , dxdt_new , t+dt );
0281
0282             //prepare dense output
0283             value_type error = prepare_dense_output( m_k_final , in , dxdt , out , dxdt_new , dt );
0284
0285             if( error > static_cast<value_type>(10) ) // we are not as accurate for interpolation as for the steps
0286             {
0287                 reject = true;
0288                 new_h = dt * pow BOOST_PREVENT_MACRO_SUBSTITUTION( error , static_cast<value_type>(-1)/(2*m_k_final+2) );
0289             } else {
0290                 t += dt;
0291             }
0292         }
0293         //set next stepsize
0294         if( !m_last_step_rejected || (new_h < dt) )
0295         {
0296             // limit step size
0297             if( m_max_dt != static_cast<time_type>(0) )
0298             {
0299                 new_h = detail::min_abs(m_max_dt, new_h);
0300             }
0301             dt = new_h;
0302         }
0303
0304         m_last_step_rejected = reject;
0305         if( reject )
0306             return fail;
0307         else
0308             return success;
0309     }
0310
0311     template< class StateType >
0312     void initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 )
0313     {
0314         m_resizer.adjust_size( x0 , detail::bind( &controlled_error_bs_type::template resize_impl< StateType > , detail::ref( *this ) , detail::_1 ) );
0315         boost::numeric::odeint::copy( x0 , get_current_state() );
0316         m_t = t0;
0317         m_dt = dt0;
0318         reset();
0319     }
0320
0321
0322     /*  =======================================================
0323      *  the actual step method that should be called from outside (maybe make try_step private?)
0324      */
0325     template< class System >
0326     std::pair< time_type , time_type > do_step( System system )
0327     {
0328         if( m_first )
0329         {
0330             typename odeint::unwrap_reference< System >::type &sys = system;
0331             sys( get_current_state() , get_current_deriv() , m_t );
0332         }
0333
0334         failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
0335         controlled_step_result res = fail;
0336         m_t_last = m_t;
0337         while( res == fail )
0338         {
0339             res = try_step( system , get_current_state() , get_current_deriv() , m_t , get_old_state() , get_old_deriv() , m_dt );
0340             m_first = false;
0341             fail_checker();  // check for overflow of failed steps
0342         }
0343         toggle_current_state();
0344         return std::make_pair( m_t_last , m_t );
0345     }
0346
0347     /* performs the interpolation from a calculated step */
0348     template< class StateOut >
0349     void calc_state( time_type t , StateOut &x ) const
0350     {
0351         do_interpolation( t , x );
0352     }
0353
0354     const state_type& current_state( void ) const
0355     {
0356         return get_current_state();
0357     }
0358
0359     time_type current_time( void ) const
0360     {
0361         return m_t;
0362     }
0363
0364     const state_type& previous_state( void ) const
0365     {
0366         return get_old_state();
0367     }
0368
0369     time_type previous_time( void ) const
0370     {
0371         return m_t_last;
0372     }
0373
0374     time_type current_time_step( void ) const
0375     {
0376         return m_dt;
0377     }
0378
0379     /** \brief Resets the internal state of the stepper. */
0380     void reset()
0381     {
0382         m_first = true;
0383         m_last_step_rejected = false;
0384     }
0385
0386     template< class StateIn >
0387     void adjust_size( const StateIn &x )
0388     {
0389         resize_impl( x );
0390         m_midpoint.adjust_size( x );
0391     }
0392
0393
0394 protected:
0395
0396     time_type m_max_dt;
0397
0398
0399 private:
0400
0401     template< class StateInOut , class StateVector >
0402     void extrapolate( size_t k , StateVector &table , const value_matrix &coeff , StateInOut &xest , size_t order_start_index = 0 )
0403     //polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
0404     {
0405         static const value_type val1( 1.0 );
0406         for( int j=k-1 ; j>0 ; --j )
0407         {
0408             m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
0409                                  typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][j + order_start_index] ,
0410                                                                                                            -coeff[k + order_start_index][j + order_start_index] ) );
0411         }
0412         m_algebra.for_each3( xest , table[0].m_v , xest ,
0413                              typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][0 + order_start_index] ,
0414                                                                                                        -coeff[k + order_start_index][0 + order_start_index]) );
0415     }
0416
0417
0418     template< class StateVector >
0419     void extrapolate_dense_out( size_t k , StateVector &table , const value_matrix &coeff , size_t order_start_index = 0 )
0420     //polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
0421     {
0422         // result is written into table[0]
0423         static const value_type val1( 1.0 );
0424         for( int j=k ; j>1 ; --j )
0425         {
0426             m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
0427                                  typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][j + order_start_index - 1] ,
0428                                                                                                            -coeff[k + order_start_index][j + order_start_index - 1] ) );
0429         }
0430         m_algebra.for_each3( table[0].m_v , table[1].m_v , table[0].m_v ,
0431                              typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][order_start_index] ,
0432                                                                                                        -coeff[k + order_start_index][order_start_index]) );
0433     }
0434
0435     time_type calc_h_opt( time_type h , value_type error , size_t k ) const
0436     {
0437         BOOST_USING_STD_MIN();
0438         BOOST_USING_STD_MAX();
0439         using std::pow;
0440
0441         value_type expo = static_cast<value_type>(1)/(m_interval_sequence[k-1]);
0442         value_type facmin = m_facmin_table[k];
0443         value_type fac;
0444         if (error == 0.0)
0445             fac = static_cast<value_type>(1)/facmin;
0446         else
0447         {
0448             fac = STEPFAC2 / pow BOOST_PREVENT_MACRO_SUBSTITUTION( error / STEPFAC1 , expo );
0449             fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( facmin/STEPFAC4 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>(static_cast<value_type>(1)/facmin) , fac ) );
0450         }
0451         return h*fac;
0452     }
0453
0454     bool in_convergence_window( size_t k ) const
0455     {
0456         if( (k == m_current_k_opt-1) && !m_last_step_rejected )
0457             return true; // decrease order only if last step was not rejected
0458         return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) );
0459     }
0460
0461     bool should_reject( value_type error , size_t k ) const
0462     {
0463         if( k == m_current_k_opt-1 )
0464         {
0465             const value_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] /
0466                 (m_interval_sequence[0]*m_interval_sequence[0]);
0467             //step will fail, criterion 17.3.17 in NR
0468             return ( error > d*d );
0469         }
0470         else if( k == m_current_k_opt )
0471         {
0472             const value_type d = m_interval_sequence[m_current_k_opt+1] / m_interval_sequence[0];
0473             return ( error > d*d );
0474         } else
0475             return error > 1.0;
0476     }
0477
0478     template< class StateIn1 , class DerivIn1 , class StateIn2 , class DerivIn2 >
0479     value_type prepare_dense_output( int k , const StateIn1 &x_start , const DerivIn1 &dxdt_start ,
0480                                      const StateIn2 & /* x_end */ , const DerivIn2 & /*dxdt_end */ , time_type dt )
0481     /* k is the order to which the result was approximated */
0482     {
0483
0484         /* compute the coefficients of the interpolation polynomial
0485          * we parametrize the interval t .. t+dt by theta = -1 .. 1
0486          * we use 2k+3 values at the interval center theta=0 to obtain the interpolation coefficients
0487          * the values are x(t+dt/2) and the derivatives dx/dt , ... d^(2k+2) x / dt^(2k+2) at the midpoints
0488          * the derivatives are approximated via finite differences
0489          * all values are obtained from interpolation of the results from the increasing orders of the midpoint calls
0490          */
0491
0492         // calculate finite difference approximations to derivatives at the midpoint
0493         for( int j = 0 ; j<=k ; j++ )
0494         {
0495             /* not working with boost units... */
0496             const value_type d = m_interval_sequence[j] / ( static_cast<value_type>(2) * dt );
0497             value_type f = 1.0; //factor 1/2 here because our interpolation interval has length 2 !!!
0498             for( int kappa = 0 ; kappa <= 2*j+1 ; ++kappa )
0499             {
0500                 calculate_finite_difference( j , kappa , f , dxdt_start );
0501                 f *= d;
0502             }
0503
0504             if( j > 0 )
0505                 extrapolate_dense_out( j , m_mp_states , m_coeff );
0506         }
0507
0508         time_type d = dt/2;
0509
0510         // extrapolate finite differences
0511         for( int kappa = 0 ; kappa<=2*k+1 ; kappa++ )
0512         {
0513             for( int j=1 ; j<=(k-kappa/2) ; ++j )
0514                 extrapolate_dense_out( j , m_diffs[kappa] , m_coeff , kappa/2 );
0515
0516             // extrapolation results are now stored in m_diffs[kappa][0]
0517
0518             // divide kappa-th derivative by kappa because we need these terms for dense output interpolation
0519             m_algebra.for_each1( m_diffs[kappa][0].m_v , typename operations_type::template scale< time_type >( static_cast<time_type>(d) ) );
0520
0521             d *= dt/(2*(kappa+2));
0522         }
0523
0524         // dense output coefficients a_0 is stored in m_mp_states[0], a_i for i = 1...2k are stored in m_diffs[i-1][0]
0525
0526         // the error is just the highest order coefficient of the interpolation polynomial
0527         // this is because we use only the midpoint theta=0 as support for the interpolation (remember that theta = -1 .. 1)
0528
0529         value_type error = 0.0;
0530         if( m_control_interpolation )
0531         {
0532             boost::numeric::odeint::copy( m_diffs[2*k+1][0].m_v , m_err.m_v );
0533             error = m_error_checker.error( m_algebra , x_start , dxdt_start , m_err.m_v , dt );
0534         }
0535
0536         return error;
0537     }
0538
0539     template< class DerivIn >
0540     void calculate_finite_difference( size_t j , size_t kappa , value_type fac , const DerivIn &dxdt )
0541     {
0542         const int m = m_interval_sequence[j]/2-1;
0543         if( kappa == 0) // no calculation required for 0th derivative of f
0544         {
0545             m_algebra.for_each2( m_diffs[0][j].m_v , m_derivs[j][m].m_v ,
0546                                  typename operations_type::template scale_sum1< value_type >( fac ) );
0547         }
0548         else
0549         {
0550             // calculate the index of m_diffs for this kappa-j-combination
0551             const int j_diffs = j - kappa/2;
0552
0553             m_algebra.for_each2( m_diffs[kappa][j_diffs].m_v , m_derivs[j][m+kappa].m_v ,
0554                                  typename operations_type::template scale_sum1< value_type >( fac ) );
0555             value_type sign = -1.0;
0556             int c = 1;
0557             //computes the j-th order finite difference for the kappa-th derivative of f at t+dt/2 using function evaluations stored in m_derivs
0558             for( int i = m+static_cast<int>(kappa)-2 ; i >= m-static_cast<int>(kappa) ; i -= 2 )
0559             {
0560                 if( i >= 0 )
0561                 {
0562                     m_algebra.for_each3( m_diffs[kappa][j_diffs].m_v , m_diffs[kappa][j_diffs].m_v , m_derivs[j][i].m_v ,
0563                                          typename operations_type::template scale_sum2< value_type , value_type >( 1.0 ,
0564                                                                                                                    sign * fac * boost::math::binomial_coefficient< value_type >( kappa , c ) ) );
0565                 }
0566                 else
0567                 {
0568                     m_algebra.for_each3( m_diffs[kappa][j_diffs].m_v , m_diffs[kappa][j_diffs].m_v , dxdt ,
0569                                          typename operations_type::template scale_sum2< value_type , value_type >( 1.0 , sign * fac ) );
0570                 }
0571                 sign *= -1;
0572                 ++c;
0573             }
0574         }
0575     }
0576
0577     template< class StateOut >
0578     void do_interpolation( time_type t , StateOut &out ) const
0579     {
0580         // interpolation polynomial is defined for theta = -1 ... 1
0581         // m_k_final is the number of order-iterations done for the last step - it governs the order of the interpolation polynomial
0582         const value_type theta = 2 * get_unit_value( (t - m_t_last) / (m_t - m_t_last) ) - 1;
0583         // we use only values at interval center, that is theta=0, for interpolation
0584         // our interpolation polynomial is thus of order 2k+2, hence we have 2k+3 terms
0585
0586         boost::numeric::odeint::copy( m_mp_states[0].m_v , out );
0587         // add remaining terms: x += a_1 theta + a2 theta^2 + ... + a_{2k} theta^{2k}
0588         value_type theta_pow( theta );
0589         for( size_t i=0 ; i<=2*m_k_final+1 ; ++i )
0590         {
0591             m_algebra.for_each3( out , out , m_diffs[i][0].m_v ,
0592                                  typename operations_type::template scale_sum2< value_type >( static_cast<value_type>(1) , theta_pow ) );
0593             theta_pow *= theta;
0594         }
0595     }
0596
0597     /* Resizer methods */
0598     template< class StateIn >
0599     bool resize_impl( const StateIn &x )
0600     {
0601         bool resized( false );
0602
0603         resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() );
0604         resized |= adjust_size_by_resizeability( m_x2 , x , typename is_resizeable<state_type>::type() );
0605         resized |= adjust_size_by_resizeability( m_dxdt1 , x , typename is_resizeable<state_type>::type() );
0606         resized |= adjust_size_by_resizeability( m_dxdt2 , x , typename is_resizeable<state_type>::type() );
0607         resized |= adjust_size_by_resizeability( m_err , x , typename is_resizeable<state_type>::type() );
0608
0609         for( size_t i = 0 ; i < m_k_max ; ++i )
0610             resized |= adjust_size_by_resizeability( m_table[i] , x , typename is_resizeable<state_type>::type() );
0611         for( size_t i = 0 ; i < m_k_max+1 ; ++i )
0612             resized |= adjust_size_by_resizeability( m_mp_states[i] , x , typename is_resizeable<state_type>::type() );
0613         for( size_t i = 0 ; i < m_k_max+1 ; ++i )
0614             for( size_t j = 0 ; j < m_derivs[i].size() ; ++j )
0615                 resized |= adjust_size_by_resizeability( m_derivs[i][j] , x , typename is_resizeable<deriv_type>::type() );
0616         for( size_t i = 0 ; i < 2*m_k_max+2 ; ++i )
0617             for( size_t j = 0 ; j < m_diffs[i].size() ; ++j )
0618                 resized |= adjust_size_by_resizeability( m_diffs[i][j] , x , typename is_resizeable<deriv_type>::type() );
0619
0620         return resized;
0621     }
0622
0623
0624     state_type& get_current_state( void )
0625     {
0626         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
0627     }
0628
0629     const state_type& get_current_state( void ) const
0630     {
0631         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
0632     }
0633
0634     state_type& get_old_state( void )
0635     {
0636         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
0637     }
0638
0639     const state_type& get_old_state( void ) const
0640     {
0641         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
0642     }
0643
0644     deriv_type& get_current_deriv( void )
0645     {
0646         return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
0647     }
0648
0649     const deriv_type& get_current_deriv( void ) const
0650     {
0651         return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
0652     }
0653
0654     deriv_type& get_old_deriv( void )
0655     {
0656         return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
0657     }
0658
0659     const deriv_type& get_old_deriv( void ) const
0660     {
0661         return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
0662     }
0663
0664
0665     void toggle_current_state( void )
0666     {
0667         m_current_state_x1 = ! m_current_state_x1;
0668     }
0669
0670
0671
0672     default_error_checker< value_type, algebra_type , operations_type > m_error_checker;
0673     modified_midpoint_dense_out< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint;
0674
0675     bool m_control_interpolation;
0676
0677     bool m_last_step_rejected;
0678     bool m_first;
0679
0680     time_type m_t;
0681     time_type m_dt;
0682     time_type m_dt_last;
0683     time_type m_t_last;
0684
0685     size_t m_current_k_opt;
0686     size_t m_k_final;
0687
0688     algebra_type m_algebra;
0689
0690     resizer_type m_resizer;
0691
0692     wrapped_state_type m_x1 , m_x2;
0693     wrapped_deriv_type m_dxdt1 , m_dxdt2;
0694     wrapped_state_type m_err;
0695     bool m_current_state_x1;
0696
0697
0698
0699     value_vector m_error; // errors of repeated midpoint steps and extrapolations
0700     int_vector m_interval_sequence; // stores the successive interval counts
0701     value_matrix m_coeff;
0702     int_vector m_cost; // costs for interval count
0703     value_vector m_facmin_table; // for precomputed facmin to save pow calls
0704
0705     state_vector_type m_table; // sequence of states for extrapolation
0706
0707     //for dense output:
0708     state_vector_type m_mp_states; // sequence of approximations of x at distance center
0709     deriv_table_type m_derivs; // table of function values
0710     deriv_table_type m_diffs; // table of function values
0711
0712     //wrapped_state_type m_a1 , m_a2 , m_a3 , m_a4;
0713
0714     value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
0715 };
0716
0717
0718
0719 /********** DOXYGEN **********/
0720
0721 /**
0722  * \class bulirsch_stoer_dense_out
0723  * \brief The Bulirsch-Stoer algorithm.
0724  *
0725  * The Bulirsch-Stoer is a controlled stepper that adjusts both step size
0726  * and order of the method. The algorithm uses the modified midpoint and
0727  * a polynomial extrapolation compute the solution. This class also provides
0728  * dense output facility.
0729  *
0730  * \tparam State The state type.
0731  * \tparam Value The value type.
0732  * \tparam Deriv The type representing the time derivative of the state.
0733  * \tparam Time The time representing the independent variable - the time.
0734  * \tparam Algebra The algebra type.
0735  * \tparam Operations The operations type.
0736  * \tparam Resizer The resizer policy type.
0737  */
0738
0739     /**
0740      * \fn bulirsch_stoer_dense_out::bulirsch_stoer_dense_out( value_type eps_abs , value_type eps_rel , value_type factor_x , value_type factor_dxdt , bool control_interpolation )
0741      * \brief Constructs the bulirsch_stoer class, including initialization of
0742      * the error bounds.
0743      *
0744      * \param eps_abs Absolute tolerance level.
0745      * \param eps_rel Relative tolerance level.
0746      * \param factor_x Factor for the weight of the state.
0747      * \param factor_dxdt Factor for the weight of the derivative.
0748      * \param control_interpolation Set true to additionally control the error of
0749      * the interpolation.
0750      */
0751
0752     /**
0753      * \fn bulirsch_stoer_dense_out::try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , DerivOut &dxdt_new , time_type &dt )
0754      * \brief Tries to perform one step.
0755      *
0756      * This method tries to do one step with step size dt. If the error estimate
0757      * is to large, the step is rejected and the method returns fail and the
0758      * step size dt is reduced. If the error estimate is acceptably small, the
0759      * step is performed, success is returned and dt might be increased to make
0760      * the steps as large as possible. This method also updates t if a step is
0761      * performed. Also, the internal order of the stepper is adjusted if required.
0762      *
0763      * \param system The system function to solve, hence the r.h.s. of the ODE.
0764      * It must fulfill the Simple System concept.
0765      * \param in The state of the ODE which should be solved.
0766      * \param dxdt The derivative of state.
0767      * \param t The value of the time. Updated if the step is successful.
0768      * \param out Used to store the result of the step.
0769      * \param dt The step size. Updated.
0770      * \return success if the step was accepted, fail otherwise.
0771      */
0772
0773     /**
0774      * \fn bulirsch_stoer_dense_out::initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 )
0775      * \brief Initializes the dense output stepper.
0776      *
0777      * \param x0 The initial state.
0778      * \param t0 The initial time.
0779      * \param dt0 The initial time step.
0780      */
0781
0782     /**
0783      * \fn bulirsch_stoer_dense_out::do_step( System system )
0784      * \brief Does one time step. This is the main method that should be used to
0785      * integrate an ODE with this stepper.
0786      * \note initialize has to be called before using this method to set the
0787      * initial conditions x,t and the stepsize.
0788      * \param system The system function to solve, hence the r.h.s. of the
0789      * ordinary differential equation. It must fulfill the Simple System concept.
0790      * \return Pair with start and end time of the integration step.
0791      */
0792
0793     /**
0794      * \fn bulirsch_stoer_dense_out::calc_state( time_type t , StateOut &x ) const
0795      * \brief Calculates the solution at an intermediate point within the last step
0796      * \param t The time at which the solution should be calculated, has to be
0797      * in the current time interval.
0798      * \param x The output variable where the result is written into.
0799      */
0800
0801     /**
0802      * \fn bulirsch_stoer_dense_out::current_state( void ) const
0803      * \brief Returns the current state of the solution.
0804      * \return The current state of the solution x(t).
0805      */
0806
0807     /**
0808      * \fn bulirsch_stoer_dense_out::current_time( void ) const
0809      * \brief Returns the current time of the solution.
0810      * \return The current time of the solution t.
0811      */
0812
0813     /**
0814      * \fn bulirsch_stoer_dense_out::previous_state( void ) const
0815      * \brief Returns the last state of the solution.
0816      * \return The last state of the solution x(t-dt).
0817      */
0818
0819     /**
0820      * \fn bulirsch_stoer_dense_out::previous_time( void ) const
0821      * \brief Returns the last time of the solution.
0822      * \return The last time of the solution t-dt.
0823      */
0824
0825     /**
0826      * \fn bulirsch_stoer_dense_out::current_time_step( void ) const
0827      * \brief Returns the current step size.
0828      * \return The current step size.
0829      */
0830
0831     /**
0832      * \fn bulirsch_stoer_dense_out::adjust_size( const StateIn &x )
0833      * \brief Adjust the size of all temporaries in the stepper manually.
0834      * \param x A state from which the size of the temporaries to be resized is deduced.
0835      */
0836
0837 }
0838 }
0839 }
0840
0841 #endif // BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED