ceph/src/boost/libs/numeric/odeint/include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp

   1 /*
   2  [auto_generated]
   3  boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
   4
   5  [begin_description]
   6  Implementaiton of the Burlish-Stoer method with dense output
   7  [end_description]
   8
   9  Copyright 2011-2015 Mario Mulansky
  10  Copyright 2011-2013 Karsten Ahnert
  11  Copyright 2012 Christoph Koke
  12
  13  Distributed under the Boost Software License, Version 1.0.
  14  (See accompanying file LICENSE_1_0.txt or
  15  copy at http://www.boost.org/LICENSE_1_0.txt)
  16  */
  17
  18
  19 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_DENSE_OUT_HPP_INCLUDED
  20 #define BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_DENSE_OUT_HPP_INCLUDED
  21
  22
  23 #include <iostream>
  24
  25 #include <algorithm>
  26
  27 #include <boost/config.hpp> // for min/max guidelines
  28
  29 #include <boost/numeric/odeint/util/bind.hpp>
  30
  31 #include <boost/math/special_functions/binomial.hpp>
  32
  33 #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
  34 #include <boost/numeric/odeint/stepper/modified_midpoint.hpp>
  35 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
  36 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
  37 #include <boost/numeric/odeint/algebra/default_operations.hpp>
  38 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
  39 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
  40
  41 #include <boost/numeric/odeint/util/state_wrapper.hpp>
  42 #include <boost/numeric/odeint/util/is_resizeable.hpp>
  43 #include <boost/numeric/odeint/util/resizer.hpp>
  44 #include <boost/numeric/odeint/util/unit_helper.hpp>
  45
  46 #include <boost/numeric/odeint/integrate/max_step_checker.hpp>
  47
  48 #include <boost/type_traits.hpp>
  49
  50
  51 namespace boost {
  52 namespace numeric {
  53 namespace odeint {
  54
  55 template<
  56     class State ,
  57     class Value = double ,
  58     class Deriv = State ,
  59     class Time = Value ,
  60     class Algebra = typename algebra_dispatcher< State >::algebra_type ,
  61     class Operations = typename operations_dispatcher< State >::operations_type ,
  62     class Resizer = initially_resizer
  63     >
  64 class bulirsch_stoer_dense_out {
  65
  66
  67 public:
  68
  69     typedef State state_type;
  70     typedef Value value_type;
  71     typedef Deriv deriv_type;
  72     typedef Time time_type;
  73     typedef Algebra algebra_type;
  74     typedef Operations operations_type;
  75     typedef Resizer resizer_type;
  76     typedef dense_output_stepper_tag stepper_category;
  77 #ifndef DOXYGEN_SKIP
  78     typedef state_wrapper< state_type > wrapped_state_type;
  79     typedef state_wrapper< deriv_type > wrapped_deriv_type;
  80
  81     typedef bulirsch_stoer_dense_out< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type;
  82
  83     typedef typename inverse_time< time_type >::type inv_time_type;
  84
  85     typedef std::vector< value_type > value_vector;
  86     typedef std::vector< time_type > time_vector;
  87     typedef std::vector< inv_time_type > inv_time_vector;  //should be 1/time_type for boost.units
  88     typedef std::vector< value_vector > value_matrix;
  89     typedef std::vector< size_t > int_vector;
  90     typedef std::vector< wrapped_state_type > state_vector_type;
  91     typedef std::vector< wrapped_deriv_type > deriv_vector_type;
  92     typedef std::vector< deriv_vector_type > deriv_table_type;
  93 #endif //DOXYGEN_SKIP
  94
  95     const static size_t m_k_max = 8;
  96
  97
  98
  99     bulirsch_stoer_dense_out(
 100         value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
 101         value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 ,
 102         time_type max_dt = static_cast<time_type>(0) ,
 103         bool control_interpolation = false )
 104         : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) ,
 105           m_max_dt(max_dt) ,
 106           m_control_interpolation( control_interpolation) ,
 107           m_last_step_rejected( false ) , m_first( true ) ,
 108           m_current_state_x1( true ) ,
 109           m_error( m_k_max ) ,
 110           m_interval_sequence( m_k_max+1 ) ,
 111           m_coeff( m_k_max+1 ) ,
 112           m_cost( m_k_max+1 ) ,
 113           m_table( m_k_max ) ,
 114           m_mp_states( m_k_max+1 ) ,
 115           m_derivs( m_k_max+1 ) ,
 116           m_diffs( 2*m_k_max+2 ) ,
 117           STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 )
 118     {
 119         BOOST_USING_STD_MIN();
 120         BOOST_USING_STD_MAX();
 121
 122         for( unsigned short i = 0; i < m_k_max+1; i++ )
 123         {
 124             /* only this specific sequence allows for dense output */
 125             m_interval_sequence[i] = 2 + 4*i;  // 2 6 10 14 ...
 126             m_derivs[i].resize( m_interval_sequence[i] );
 127             if( i == 0 )
 128                 m_cost[i] = m_interval_sequence[i];
 129             else
 130                 m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
 131             m_coeff[i].resize(i);
 132             for( size_t k = 0 ; k < i ; ++k  )
 133             {
 134                 const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
 135                 m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
 136             }
 137             // crude estimate of optimal order
 138
 139             m_current_k_opt = 4;
 140             /* no calculation because log10 might not exist for value_type!
 141             const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >( 1.0E-12 ) ) ) * 0.6 + 0.5 );
 142             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 ) ));
 143             */
 144         }
 145         int num = 1;
 146         for( int i = 2*(m_k_max)+1 ; i >=0  ; i-- )
 147         {
 148             m_diffs[i].resize( num );
 149             num += (i+1)%2;
 150         }
 151     }
 152
 153     template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
 154     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , DerivOut &dxdt_new , time_type &dt )
 155     {
 156         if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
 157         {
 158             // given step size is bigger then max_dt
 159             // set limit and return fail
 160             dt = m_max_dt;
 161             return fail;
 162         }
 163
 164         BOOST_USING_STD_MIN();
 165         BOOST_USING_STD_MAX();
 166         using std::pow;
 167
 168         static const value_type val1( 1.0 );
 169
 170         bool reject( true );
 171
 172         time_vector h_opt( m_k_max+1 );
 173         inv_time_vector work( m_k_max+1 );
 174
 175         m_k_final = 0;
 176         time_type new_h = dt;
 177
 178         //std::cout << "t=" << t <<", dt=" << dt << ", k_opt=" << m_current_k_opt << ", first: " << m_first << std::endl;
 179
 180         for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ )
 181         {
 182             m_midpoint.set_steps( m_interval_sequence[k] );
 183             if( k == 0 )
 184             {
 185                 m_midpoint.do_step( system , in , dxdt , t , out , dt , m_mp_states[k].m_v , m_derivs[k]);
 186             }
 187             else
 188             {
 189                 m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt , m_mp_states[k].m_v , m_derivs[k] );
 190                 extrapolate( k , m_table , m_coeff , out );
 191                 // get error estimate
 192                 m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v ,
 193                                      typename operations_type::template scale_sum2< value_type , value_type >( val1 , -val1 ) );
 194                 const value_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt );
 195                 h_opt[k] = calc_h_opt( dt , error , k );
 196                 work[k] = static_cast<value_type>( m_cost[k] ) / h_opt[k];
 197
 198                 m_k_final = k;
 199
 200                 if( (k == m_current_k_opt-1) || m_first )
 201                 { // convergence before k_opt ?
 202                     if( error < 1.0 )
 203                     {
 204                         //convergence
 205                         reject = false;
 206                         if( (work[k] < KFAC2*work[k-1]) || (m_current_k_opt <= 2) )
 207                         {
 208                             // leave order as is (except we were in first round)
 209                             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 ) );
 210                             new_h = h_opt[k] * static_cast<value_type>( m_cost[k+1] ) / static_cast<value_type>( m_cost[k] );
 211                         } else {
 212                             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) ) );
 213                             new_h = h_opt[k];
 214                         }
 215                         break;
 216                     }
 217                     else if( should_reject( error , k ) && !m_first )
 218                     {
 219                         reject = true;
 220                         new_h = h_opt[k];
 221                         break;
 222                     }
 223                 }
 224                 if( k == m_current_k_opt )
 225                 { // convergence at k_opt ?
 226                     if( error < 1.0 )
 227                     {
 228                         //convergence
 229                         reject = false;
 230                         if( (work[k-1] < KFAC2*work[k]) )
 231                         {
 232                             m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
 233                             new_h = h_opt[m_current_k_opt];
 234                         }
 235                         else if( (work[k] < KFAC2*work[k-1]) && !m_last_step_rejected )
 236                         {
 237                             m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , static_cast<int>(m_current_k_opt)+1 );
 238                             new_h = h_opt[k]*static_cast<value_type>( m_cost[m_current_k_opt] ) / static_cast<value_type>( m_cost[k] );
 239                         } else
 240                             new_h = h_opt[m_current_k_opt];
 241                         break;
 242                     }
 243                     else if( should_reject( error , k ) )
 244                     {
 245                         reject = true;
 246                         new_h = h_opt[m_current_k_opt];
 247                         break;
 248                     }
 249                 }
 250                 if( k == m_current_k_opt+1 )
 251                 { // convergence at k_opt+1 ?
 252                     if( error < 1.0 )
 253                     {   //convergence
 254                         reject = false;
 255                         if( work[k-2] < KFAC2*work[k-1] )
 256                             m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
 257                         if( (work[k] < KFAC2*work[m_current_k_opt]) && !m_last_step_rejected )
 258                             m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , static_cast<int>(k) );
 259                         new_h = h_opt[m_current_k_opt];
 260                     } else
 261                     {
 262                         reject = true;
 263                         new_h = h_opt[m_current_k_opt];
 264                     }
 265                     break;
 266                 }
 267             }
 268         }
 269
 270         if( !reject )
 271         {
 272
 273             //calculate dxdt for next step and dense output
 274             typename odeint::unwrap_reference< System >::type &sys = system;
 275             sys( out , dxdt_new , t+dt );
 276
 277             //prepare dense output
 278             value_type error = prepare_dense_output( m_k_final , in , dxdt , out , dxdt_new , dt );
 279
 280             if( error > static_cast<value_type>(10) ) // we are not as accurate for interpolation as for the steps
 281             {
 282                 reject = true;
 283                 new_h = dt * pow BOOST_PREVENT_MACRO_SUBSTITUTION( error , static_cast<value_type>(-1)/(2*m_k_final+2) );
 284             } else {
 285                 t += dt;
 286             }
 287         }
 288         //set next stepsize
 289         if( !m_last_step_rejected || (new_h < dt) )
 290         {
 291             // limit step size
 292             if( m_max_dt != static_cast<time_type>(0) )
 293             {
 294                 new_h = detail::min_abs(m_max_dt, new_h);
 295             }
 296             dt = new_h;
 297         }
 298
 299         m_last_step_rejected = reject;
 300         if( reject )
 301             return fail;
 302         else
 303             return success;
 304     }
 305
 306     template< class StateType >
 307     void initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 )
 308     {
 309         m_resizer.adjust_size( x0 , detail::bind( &controlled_error_bs_type::template resize_impl< StateType > , detail::ref( *this ) , detail::_1 ) );
 310         boost::numeric::odeint::copy( x0 , get_current_state() );
 311         m_t = t0;
 312         m_dt = dt0;
 313         reset();
 314     }
 315
 316
 317     /*  =======================================================
 318      *  the actual step method that should be called from outside (maybe make try_step private?)
 319      */
 320     template< class System >
 321     std::pair< time_type , time_type > do_step( System system )
 322     {
 323         if( m_first )
 324         {
 325             typename odeint::unwrap_reference< System >::type &sys = system;
 326             sys( get_current_state() , get_current_deriv() , m_t );
 327         }
 328
 329         failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
 330         controlled_step_result res = fail;
 331         m_t_last = m_t;
 332         while( res == fail )
 333         {
 334             res = try_step( system , get_current_state() , get_current_deriv() , m_t , get_old_state() , get_old_deriv() , m_dt );
 335             m_first = false;
 336             fail_checker();  // check for overflow of failed steps
 337         }
 338         toggle_current_state();
 339         return std::make_pair( m_t_last , m_t );
 340     }
 341
 342     /* performs the interpolation from a calculated step */
 343     template< class StateOut >
 344     void calc_state( time_type t , StateOut &x ) const
 345     {
 346         do_interpolation( t , x );
 347     }
 348
 349     const state_type& current_state( void ) const
 350     {
 351         return get_current_state();
 352     }
 353
 354     time_type current_time( void ) const
 355     {
 356         return m_t;
 357     }
 358
 359     const state_type& previous_state( void ) const
 360     {
 361         return get_old_state();
 362     }
 363
 364     time_type previous_time( void ) const
 365     {
 366         return m_t_last;
 367     }
 368
 369     time_type current_time_step( void ) const
 370     {
 371         return m_dt;
 372     }
 373
 374     /** \brief Resets the internal state of the stepper. */
 375     void reset()
 376     {
 377         m_first = true;
 378         m_last_step_rejected = false;
 379     }
 380
 381     template< class StateIn >
 382     void adjust_size( const StateIn &x )
 383     {
 384         resize_impl( x );
 385         m_midpoint.adjust_size( x );
 386     }
 387
 388
 389 private:
 390
 391     template< class StateInOut , class StateVector >
 392     void extrapolate( size_t k , StateVector &table , const value_matrix &coeff , StateInOut &xest , size_t order_start_index = 0 )
 393     //polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
 394     {
 395         static const value_type val1( 1.0 );
 396         for( int j=k-1 ; j>0 ; --j )
 397         {
 398             m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
 399                                  typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][j + order_start_index] ,
 400                                                                                                            -coeff[k + order_start_index][j + order_start_index] ) );
 401         }
 402         m_algebra.for_each3( xest , table[0].m_v , xest ,
 403                              typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][0 + order_start_index] ,
 404                                                                                                        -coeff[k + order_start_index][0 + order_start_index]) );
 405     }
 406
 407
 408     template< class StateVector >
 409     void extrapolate_dense_out( size_t k , StateVector &table , const value_matrix &coeff , size_t order_start_index = 0 )
 410     //polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
 411     {
 412         // result is written into table[0]
 413         static const value_type val1( 1.0 );
 414         for( int j=k ; j>1 ; --j )
 415         {
 416             m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
 417                                  typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][j + order_start_index - 1] ,
 418                                                                                                            -coeff[k + order_start_index][j + order_start_index - 1] ) );
 419         }
 420         m_algebra.for_each3( table[0].m_v , table[1].m_v , table[0].m_v ,
 421                              typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][order_start_index] ,
 422                                                                                                        -coeff[k + order_start_index][order_start_index]) );
 423     }
 424
 425     time_type calc_h_opt( time_type h , value_type error , size_t k ) const
 426     {
 427         BOOST_USING_STD_MIN();
 428         BOOST_USING_STD_MAX();
 429         using std::pow;
 430
 431         value_type expo = static_cast<value_type>(1)/(m_interval_sequence[k-1]);
 432         value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
 433         value_type fac;
 434         if (error == 0.0)
 435             fac = static_cast<value_type>(1)/facmin;
 436         else
 437         {
 438             fac = STEPFAC2 / pow BOOST_PREVENT_MACRO_SUBSTITUTION( error / STEPFAC1 , expo );
 439             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 ) );
 440         }
 441         return h*fac;
 442     }
 443
 444     bool in_convergence_window( size_t k ) const
 445     {
 446         if( (k == m_current_k_opt-1) && !m_last_step_rejected )
 447             return true; // decrease order only if last step was not rejected
 448         return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) );
 449     }
 450
 451     bool should_reject( value_type error , size_t k ) const
 452     {
 453         if( k == m_current_k_opt-1 )
 454         {
 455             const value_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] /
 456                 (m_interval_sequence[0]*m_interval_sequence[0]);
 457             //step will fail, criterion 17.3.17 in NR
 458             return ( error > d*d );
 459         }
 460         else if( k == m_current_k_opt )
 461         {
 462             const value_type d = m_interval_sequence[m_current_k_opt+1] / m_interval_sequence[0];
 463             return ( error > d*d );
 464         } else
 465             return error > 1.0;
 466     }
 467
 468     template< class StateIn1 , class DerivIn1 , class StateIn2 , class DerivIn2 >
 469     value_type prepare_dense_output( int k , const StateIn1 &x_start , const DerivIn1 &dxdt_start ,
 470                                      const StateIn2 & /* x_end */ , const DerivIn2 & /*dxdt_end */ , time_type dt )
 471     /* k is the order to which the result was approximated */
 472     {
 473
 474         /* compute the coefficients of the interpolation polynomial
 475          * we parametrize the interval t .. t+dt by theta = -1 .. 1
 476          * we use 2k+3 values at the interval center theta=0 to obtain the interpolation coefficients
 477          * the values are x(t+dt/2) and the derivatives dx/dt , ... d^(2k+2) x / dt^(2k+2) at the midpoints
 478          * the derivatives are approximated via finite differences
 479          * all values are obtained from interpolation of the results from the increasing orders of the midpoint calls
 480          */
 481
 482         // calculate finite difference approximations to derivatives at the midpoint
 483         for( int j = 0 ; j<=k ; j++ )
 484         {
 485             /* not working with boost units... */
 486             const value_type d = m_interval_sequence[j] / ( static_cast<value_type>(2) * dt );
 487             value_type f = 1.0; //factor 1/2 here because our interpolation interval has length 2 !!!
 488             for( int kappa = 0 ; kappa <= 2*j+1 ; ++kappa )
 489             {
 490                 calculate_finite_difference( j , kappa , f , dxdt_start );
 491                 f *= d;
 492             }
 493
 494             if( j > 0 )
 495                 extrapolate_dense_out( j , m_mp_states , m_coeff );
 496         }
 497
 498         time_type d = dt/2;
 499
 500         // extrapolate finite differences
 501         for( int kappa = 0 ; kappa<=2*k+1 ; kappa++ )
 502         {
 503             for( int j=1 ; j<=(k-kappa/2) ; ++j )
 504                 extrapolate_dense_out( j , m_diffs[kappa] , m_coeff , kappa/2 );
 505
 506             // extrapolation results are now stored in m_diffs[kappa][0]
 507
 508             // divide kappa-th derivative by kappa because we need these terms for dense output interpolation
 509             m_algebra.for_each1( m_diffs[kappa][0].m_v , typename operations_type::template scale< time_type >( static_cast<time_type>(d) ) );
 510
 511             d *= dt/(2*(kappa+2));
 512         }
 513
 514         // 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]
 515
 516         // the error is just the highest order coefficient of the interpolation polynomial
 517         // this is because we use only the midpoint theta=0 as support for the interpolation (remember that theta = -1 .. 1)
 518
 519         value_type error = 0.0;
 520         if( m_control_interpolation )
 521         {
 522             boost::numeric::odeint::copy( m_diffs[2*k+1][0].m_v , m_err.m_v );
 523             error = m_error_checker.error( m_algebra , x_start , dxdt_start , m_err.m_v , dt );
 524         }
 525
 526         return error;
 527     }
 528
 529     template< class DerivIn >
 530     void calculate_finite_difference( size_t j , size_t kappa , value_type fac , const DerivIn &dxdt )
 531     {
 532         const int m = m_interval_sequence[j]/2-1;
 533         if( kappa == 0) // no calculation required for 0th derivative of f
 534         {
 535             m_algebra.for_each2( m_diffs[0][j].m_v , m_derivs[j][m].m_v ,
 536                                  typename operations_type::template scale_sum1< value_type >( fac ) );
 537         }
 538         else
 539         {
 540             // calculate the index of m_diffs for this kappa-j-combination
 541             const int j_diffs = j - kappa/2;
 542
 543             m_algebra.for_each2( m_diffs[kappa][j_diffs].m_v , m_derivs[j][m+kappa].m_v ,
 544                                  typename operations_type::template scale_sum1< value_type >( fac ) );
 545             value_type sign = -1.0;
 546             int c = 1;
 547             //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
 548             for( int i = m+static_cast<int>(kappa)-2 ; i >= m-static_cast<int>(kappa) ; i -= 2 )
 549             {
 550                 if( i >= 0 )
 551                 {
 552                     m_algebra.for_each3( m_diffs[kappa][j_diffs].m_v , m_diffs[kappa][j_diffs].m_v , m_derivs[j][i].m_v ,
 553                                          typename operations_type::template scale_sum2< value_type , value_type >( 1.0 ,
 554                                                                                                                    sign * fac * boost::math::binomial_coefficient< value_type >( kappa , c ) ) );
 555                 }
 556                 else
 557                 {
 558                     m_algebra.for_each3( m_diffs[kappa][j_diffs].m_v , m_diffs[kappa][j_diffs].m_v , dxdt ,
 559                                          typename operations_type::template scale_sum2< value_type , value_type >( 1.0 , sign * fac ) );
 560                 }
 561                 sign *= -1;
 562                 ++c;
 563             }
 564         }
 565     }
 566
 567     template< class StateOut >
 568     void do_interpolation( time_type t , StateOut &out ) const
 569     {
 570         // interpolation polynomial is defined for theta = -1 ... 1
 571         // m_k_final is the number of order-iterations done for the last step - it governs the order of the interpolation polynomial
 572         const value_type theta = 2 * get_unit_value( (t - m_t_last) / (m_t - m_t_last) ) - 1;
 573         // we use only values at interval center, that is theta=0, for interpolation
 574         // our interpolation polynomial is thus of order 2k+2, hence we have 2k+3 terms
 575
 576         boost::numeric::odeint::copy( m_mp_states[0].m_v , out );
 577         // add remaining terms: x += a_1 theta + a2 theta^2 + ... + a_{2k} theta^{2k}
 578         value_type theta_pow( theta );
 579         for( size_t i=0 ; i<=2*m_k_final+1 ; ++i )
 580         {
 581             m_algebra.for_each3( out , out , m_diffs[i][0].m_v ,
 582                                  typename operations_type::template scale_sum2< value_type >( static_cast<value_type>(1) , theta_pow ) );
 583             theta_pow *= theta;
 584         }
 585     }
 586
 587     /* Resizer methods */
 588     template< class StateIn >
 589     bool resize_impl( const StateIn &x )
 590     {
 591         bool resized( false );
 592
 593         resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() );
 594         resized |= adjust_size_by_resizeability( m_x2 , x , typename is_resizeable<state_type>::type() );
 595         resized |= adjust_size_by_resizeability( m_dxdt1 , x , typename is_resizeable<state_type>::type() );
 596         resized |= adjust_size_by_resizeability( m_dxdt2 , x , typename is_resizeable<state_type>::type() );
 597         resized |= adjust_size_by_resizeability( m_err , x , typename is_resizeable<state_type>::type() );
 598
 599         for( size_t i = 0 ; i < m_k_max ; ++i )
 600             resized |= adjust_size_by_resizeability( m_table[i] , x , typename is_resizeable<state_type>::type() );
 601         for( size_t i = 0 ; i < m_k_max+1 ; ++i )
 602             resized |= adjust_size_by_resizeability( m_mp_states[i] , x , typename is_resizeable<state_type>::type() );
 603         for( size_t i = 0 ; i < m_k_max+1 ; ++i )
 604             for( size_t j = 0 ; j < m_derivs[i].size() ; ++j )
 605                 resized |= adjust_size_by_resizeability( m_derivs[i][j] , x , typename is_resizeable<deriv_type>::type() );
 606         for( size_t i = 0 ; i < 2*m_k_max+2 ; ++i )
 607             for( size_t j = 0 ; j < m_diffs[i].size() ; ++j )
 608                 resized |= adjust_size_by_resizeability( m_diffs[i][j] , x , typename is_resizeable<deriv_type>::type() );
 609
 610         return resized;
 611     }
 612
 613
 614     state_type& get_current_state( void )
 615     {
 616         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
 617     }
 618
 619     const state_type& get_current_state( void ) const
 620     {
 621         return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
 622     }
 623
 624     state_type& get_old_state( void )
 625     {
 626         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
 627     }
 628
 629     const state_type& get_old_state( void ) const
 630     {
 631         return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
 632     }
 633
 634     deriv_type& get_current_deriv( void )
 635     {
 636         return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
 637     }
 638
 639     const deriv_type& get_current_deriv( void ) const
 640     {
 641         return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
 642     }
 643
 644     deriv_type& get_old_deriv( void )
 645     {
 646         return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
 647     }
 648
 649     const deriv_type& get_old_deriv( void ) const
 650     {
 651         return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
 652     }
 653
 654
 655     void toggle_current_state( void )
 656     {
 657         m_current_state_x1 = ! m_current_state_x1;
 658     }
 659
 660
 661
 662     default_error_checker< value_type, algebra_type , operations_type > m_error_checker;
 663     modified_midpoint_dense_out< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint;
 664
 665     time_type m_max_dt;
 666
 667     bool m_control_interpolation;
 668
 669     bool m_last_step_rejected;
 670     bool m_first;
 671
 672     time_type m_t;
 673     time_type m_dt;
 674     time_type m_dt_last;
 675     time_type m_t_last;
 676
 677     size_t m_current_k_opt;
 678     size_t m_k_final;
 679
 680     algebra_type m_algebra;
 681
 682     resizer_type m_resizer;
 683
 684     wrapped_state_type m_x1 , m_x2;
 685     wrapped_deriv_type m_dxdt1 , m_dxdt2;
 686     wrapped_state_type m_err;
 687     bool m_current_state_x1;
 688
 689
 690
 691     value_vector m_error; // errors of repeated midpoint steps and extrapolations
 692     int_vector m_interval_sequence; // stores the successive interval counts
 693     value_matrix m_coeff;
 694     int_vector m_cost; // costs for interval count
 695
 696     state_vector_type m_table; // sequence of states for extrapolation
 697
 698     //for dense output:
 699     state_vector_type m_mp_states; // sequence of approximations of x at distance center
 700     deriv_table_type m_derivs; // table of function values
 701     deriv_table_type m_diffs; // table of function values
 702
 703     //wrapped_state_type m_a1 , m_a2 , m_a3 , m_a4;
 704
 705     value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
 706 };
 707
 708
 709
 710 /********** DOXYGEN **********/
 711
 712 /**
 713  * \class bulirsch_stoer_dense_out
 714  * \brief The Bulirsch-Stoer algorithm.
 715  *
 716  * The Bulirsch-Stoer is a controlled stepper that adjusts both step size
 717  * and order of the method. The algorithm uses the modified midpoint and
 718  * a polynomial extrapolation compute the solution. This class also provides
 719  * dense output facility.
 720  *
 721  * \tparam State The state type.
 722  * \tparam Value The value type.
 723  * \tparam Deriv The type representing the time derivative of the state.
 724  * \tparam Time The time representing the independent variable - the time.
 725  * \tparam Algebra The algebra type.
 726  * \tparam Operations The operations type.
 727  * \tparam Resizer The resizer policy type.
 728  */
 729
 730     /**
 731      * \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 )
 732      * \brief Constructs the bulirsch_stoer class, including initialization of
 733      * the error bounds.
 734      *
 735      * \param eps_abs Absolute tolerance level.
 736      * \param eps_rel Relative tolerance level.
 737      * \param factor_x Factor for the weight of the state.
 738      * \param factor_dxdt Factor for the weight of the derivative.
 739      * \param control_interpolation Set true to additionally control the error of
 740      * the interpolation.
 741      */
 742
 743     /**
 744      * \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 )
 745      * \brief Tries to perform one step.
 746      *
 747      * This method tries to do one step with step size dt. If the error estimate
 748      * is to large, the step is rejected and the method returns fail and the
 749      * step size dt is reduced. If the error estimate is acceptably small, the
 750      * step is performed, success is returned and dt might be increased to make
 751      * the steps as large as possible. This method also updates t if a step is
 752      * performed. Also, the internal order of the stepper is adjusted if required.
 753      *
 754      * \param system The system function to solve, hence the r.h.s. of the ODE.
 755      * It must fulfill the Simple System concept.
 756      * \param in The state of the ODE which should be solved.
 757      * \param dxdt The derivative of state.
 758      * \param t The value of the time. Updated if the step is successful.
 759      * \param out Used to store the result of the step.
 760      * \param dt The step size. Updated.
 761      * \return success if the step was accepted, fail otherwise.
 762      */
 763
 764     /**
 765      * \fn bulirsch_stoer_dense_out::initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 )
 766      * \brief Initializes the dense output stepper.
 767      *
 768      * \param x0 The initial state.
 769      * \param t0 The initial time.
 770      * \param dt0 The initial time step.
 771      */
 772
 773     /**
 774      * \fn bulirsch_stoer_dense_out::do_step( System system )
 775      * \brief Does one time step. This is the main method that should be used to
 776      * integrate an ODE with this stepper.
 777      * \note initialize has to be called before using this method to set the
 778      * initial conditions x,t and the stepsize.
 779      * \param system The system function to solve, hence the r.h.s. of the
 780      * ordinary differential equation. It must fulfill the Simple System concept.
 781      * \return Pair with start and end time of the integration step.
 782      */
 783
 784     /**
 785      * \fn bulirsch_stoer_dense_out::calc_state( time_type t , StateOut &x ) const
 786      * \brief Calculates the solution at an intermediate point within the last step
 787      * \param t The time at which the solution should be calculated, has to be
 788      * in the current time interval.
 789      * \param x The output variable where the result is written into.
 790      */
 791
 792     /**
 793      * \fn bulirsch_stoer_dense_out::current_state( void ) const
 794      * \brief Returns the current state of the solution.
 795      * \return The current state of the solution x(t).
 796      */
 797
 798     /**
 799      * \fn bulirsch_stoer_dense_out::current_time( void ) const
 800      * \brief Returns the current time of the solution.
 801      * \return The current time of the solution t.
 802      */
 803
 804     /**
 805      * \fn bulirsch_stoer_dense_out::previous_state( void ) const
 806      * \brief Returns the last state of the solution.
 807      * \return The last state of the solution x(t-dt).
 808      */
 809
 810     /**
 811      * \fn bulirsch_stoer_dense_out::previous_time( void ) const
 812      * \brief Returns the last time of the solution.
 813      * \return The last time of the solution t-dt.
 814      */
 815
 816     /**
 817      * \fn bulirsch_stoer_dense_out::current_time_step( void ) const
 818      * \brief Returns the current step size.
 819      * \return The current step size.
 820      */
 821
 822     /**
 823      * \fn bulirsch_stoer_dense_out::adjust_size( const StateIn &x )
 824      * \brief Adjust the size of all temporaries in the stepper manually.
 825      * \param x A state from which the size of the temporaries to be resized is deduced.
 826      */
 827
 828 }
 829 }
 830 }
 831
 832 #endif // BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED