ceph/src/boost/boost/numeric/odeint/stepper/bulirsch_stoer.hpp

   1 /*
   2   [auto_generated]
   3   boost/numeric/odeint/stepper/bulirsch_stoer.hpp
   4
   5   [begin_description]
   6   Implementation of the Burlish-Stoer method. As described in
   7   Ernst Hairer, Syvert Paul Norsett, Gerhard Wanner
   8   Solving Ordinary Differential Equations I. Nonstiff Problems.
   9   Springer Series in Comput. Mathematics, Vol. 8, Springer-Verlag 1987, Second revised edition 1993.
  10   [end_description]
  11
  12   Copyright 2011-2013 Mario Mulansky
  13   Copyright 2011-2013 Karsten Ahnert
  14   Copyright 2012 Christoph Koke
  15
  16   Distributed under the Boost Software License, Version 1.0.
  17   (See accompanying file LICENSE_1_0.txt or
  18   copy at http://www.boost.org/LICENSE_1_0.txt)
  19 */
  20
  21
  22 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED
  23 #define BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED
  24
  25
  26 #include <iostream>
  27
  28 #include <algorithm>
  29
  30 #include <boost/config.hpp> // for min/max guidelines
  31
  32 #include <boost/numeric/odeint/util/bind.hpp>
  33 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
  34
  35 #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
  36 #include <boost/numeric/odeint/stepper/modified_midpoint.hpp>
  37 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
  38 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
  39 #include <boost/numeric/odeint/algebra/default_operations.hpp>
  40 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
  41 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
  42
  43 #include <boost/numeric/odeint/util/state_wrapper.hpp>
  44 #include <boost/numeric/odeint/util/is_resizeable.hpp>
  45 #include <boost/numeric/odeint/util/resizer.hpp>
  46 #include <boost/numeric/odeint/util/unit_helper.hpp>
  47 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
  48
  49 namespace boost {
  50 namespace numeric {
  51 namespace odeint {
  52
  53 template<
  54     class State ,
  55     class Value = double ,
  56     class Deriv = State ,
  57     class Time = Value ,
  58     class Algebra = typename algebra_dispatcher< State >::algebra_type ,
  59     class Operations = typename operations_dispatcher< State >::operations_type ,
  60     class Resizer = initially_resizer
  61     >
  62 class bulirsch_stoer {
  63
  64 public:
  65
  66     typedef State state_type;
  67     typedef Value value_type;
  68     typedef Deriv deriv_type;
  69     typedef Time time_type;
  70     typedef Algebra algebra_type;
  71     typedef Operations operations_type;
  72     typedef Resizer resizer_type;
  73 #ifndef DOXYGEN_SKIP
  74     typedef state_wrapper< state_type > wrapped_state_type;
  75     typedef state_wrapper< deriv_type > wrapped_deriv_type;
  76     typedef controlled_stepper_tag stepper_category;
  77
  78     typedef bulirsch_stoer< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type;
  79
  80     typedef typename inverse_time< time_type >::type inv_time_type;
  81
  82     typedef std::vector< value_type > value_vector;
  83     typedef std::vector< time_type > time_vector;
  84     typedef std::vector< inv_time_type > inv_time_vector;  //should be 1/time_type for boost.units
  85     typedef std::vector< value_vector > value_matrix;
  86     typedef std::vector< size_t > int_vector;
  87     typedef std::vector< wrapped_state_type > state_table_type;
  88 #endif //DOXYGEN_SKIP
  89     const static size_t m_k_max = 8;
  90
  91     bulirsch_stoer(
  92         value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
  93         value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 ,
  94         time_type max_dt = static_cast<time_type>(0))
  95         : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , m_midpoint() ,
  96           m_last_step_rejected( false ) , m_first( true ) ,
  97           m_max_dt(max_dt) ,
  98           m_interval_sequence( m_k_max+1 ) ,
  99           m_coeff( m_k_max+1 ) ,
 100           m_cost( m_k_max+1 ) ,
 101           m_facmin_table( m_k_max+1 ) ,
 102           m_table( m_k_max ) ,
 103           STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 )
 104     {
 105         BOOST_USING_STD_MIN();
 106         BOOST_USING_STD_MAX();
 107         /* initialize sequence of stage numbers and work */
 108         for( unsigned short i = 0; i < m_k_max+1; i++ )
 109         {
 110             m_interval_sequence[i] = 2 * (i+1);
 111             if( i == 0 )
 112                 m_cost[i] = m_interval_sequence[i];
 113             else
 114                 m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
 115             m_coeff[i].resize(i);
 116             m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) );
 117             for( size_t k = 0 ; k < i ; ++k  )
 118             {
 119                 const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
 120                 m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
 121             }
 122         }
 123         reset();
 124     }
 125
 126
 127     /*
 128      * Version 1 : try_step( sys , x , t , dt )
 129      *
 130      * The overloads are needed to solve the forwarding problem
 131      */
 132     template< class System , class StateInOut >
 133     controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
 134     {
 135         return try_step_v1( system , x , t, dt );
 136     }
 137
 138     /**
 139      * \brief Second version to solve the forwarding problem, can be used with Boost.Range as StateInOut.
 140      */
 141     template< class System , class StateInOut >
 142     controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
 143     {
 144         return try_step_v1( system , x , t, dt );
 145     }
 146
 147     /*
 148      * Version 2 : try_step( sys , x , dxdt , t , dt )
 149      *
 150      * this version does not solve the forwarding problem, boost.range can not be used
 151      */
 152     template< class System , class StateInOut , class DerivIn >
 153     controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt )
 154     {
 155         m_xnew_resizer.adjust_size( x , detail::bind( &controlled_error_bs_type::template resize_m_xnew< StateInOut > , detail::ref( *this ) , detail::_1 ) );
 156         controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt );
 157         if( res == success )
 158         {
 159             boost::numeric::odeint::copy( m_xnew.m_v , x );
 160         }
 161         return res;
 162     }
 163
 164     /*
 165      * Version 3 : try_step( sys , in , t , out , dt )
 166      *
 167      * this version does not solve the forwarding problem, boost.range can not be used
 168      */
 169     template< class System , class StateIn , class StateOut >
 170     typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type
 171     try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
 172     {
 173         typename odeint::unwrap_reference< System >::type &sys = system;
 174         m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_m_dxdt< StateIn > , detail::ref( *this ) , detail::_1 ) );
 175         sys( in , m_dxdt.m_v , t );
 176         return try_step( system , in , m_dxdt.m_v , t , out , dt );
 177     }
 178
 179
 180     /*
 181      * Full version : try_step( sys , in , dxdt_in , t , out , dt )
 182      *
 183      * contains the actual implementation
 184      */
 185     template< class System , class StateIn , class DerivIn , class StateOut >
 186     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
 187     {
 188         if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
 189         {
 190             // given step size is bigger then max_dt
 191             // set limit and return fail
 192             dt = m_max_dt;
 193             return fail;
 194         }
 195
 196         BOOST_USING_STD_MIN();
 197         BOOST_USING_STD_MAX();
 198
 199         static const value_type val1( 1.0 );
 200
 201         if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) )
 202         {
 203             reset(); // system resized -> reset
 204         }
 205
 206         if( dt != m_dt_last )
 207         {
 208             reset(); // step size changed from outside -> reset
 209         }
 210
 211         bool reject( true );
 212
 213         time_vector h_opt( m_k_max+1 );
 214         inv_time_vector work( m_k_max+1 );
 215
 216         time_type new_h = dt;
 217
 218         /* m_current_k_opt is the estimated current optimal stage number */
 219         for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ )
 220         {
 221             /* the stage counts are stored in m_interval_sequence */
 222             m_midpoint.set_steps( m_interval_sequence[k] );
 223             if( k == 0 )
 224             {
 225                 m_midpoint.do_step( system , in , dxdt , t , out , dt );
 226                 /* the first step, nothing more to do */
 227             }
 228             else
 229             {
 230                 m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt );
 231                 extrapolate( k , m_table , m_coeff , out );
 232                 // get error estimate
 233                 m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v ,
 234                                      typename operations_type::template scale_sum2< value_type , value_type >( val1 , -val1 ) );
 235                 const value_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt );
 236                 h_opt[k] = calc_h_opt( dt , error , k );
 237                 work[k] = static_cast<value_type>( m_cost[k] ) / h_opt[k];
 238
 239                 if( (k == m_current_k_opt-1) || m_first )
 240                 { // convergence before k_opt ?
 241                     if( error < 1.0 )
 242                     {
 243                         //convergence
 244                         reject = false;
 245                         if( (work[k] < KFAC2*work[k-1]) || (m_current_k_opt <= 2) )
 246                         {
 247                             // leave order as is (except we were in first round)
 248                             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 ) );
 249                             new_h = h_opt[k];
 250                             new_h *= static_cast<value_type>( m_cost[k+1] ) / static_cast<value_type>( m_cost[k] );
 251                         } else {
 252                             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) ) );
 253                             new_h = h_opt[k];
 254                         }
 255                         break;
 256                     }
 257                     else if( should_reject( error , k ) && !m_first )
 258                     {
 259                         reject = true;
 260                         new_h = h_opt[k];
 261                         break;
 262                     }
 263                 }
 264                 if( k == m_current_k_opt )
 265                 { // convergence at k_opt ?
 266                     if( error < 1.0 )
 267                     {
 268                         //convergence
 269                         reject = false;
 270                         if( (work[k-1] < KFAC2*work[k]) )
 271                         {
 272                             m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
 273                             new_h = h_opt[m_current_k_opt];
 274                         }
 275                         else if( (work[k] < KFAC2*work[k-1]) && !m_last_step_rejected )
 276                         {
 277                             m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max-1) , static_cast<int>(m_current_k_opt)+1 );
 278                             new_h = h_opt[k];
 279                             new_h *= static_cast<value_type>(m_cost[m_current_k_opt])/static_cast<value_type>(m_cost[k]);
 280                         } else
 281                             new_h = h_opt[m_current_k_opt];
 282                         break;
 283                     }
 284                     else if( should_reject( error , k ) )
 285                     {
 286                         reject = true;
 287                         new_h = h_opt[m_current_k_opt];
 288                         break;
 289                     }
 290                 }
 291                 if( k == m_current_k_opt+1 )
 292                 { // convergence at k_opt+1 ?
 293                     if( error < 1.0 )
 294                     {   //convergence
 295                         reject = false;
 296                         if( work[k-2] < KFAC2*work[k-1] )
 297                             m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
 298                         if( (work[k] < KFAC2*work[m_current_k_opt]) && !m_last_step_rejected )
 299                             m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , static_cast<int>(k) );
 300                         new_h = h_opt[m_current_k_opt];
 301                     } else
 302                     {
 303                         reject = true;
 304                         new_h = h_opt[m_current_k_opt];
 305                     }
 306                     break;
 307                 }
 308             }
 309         }
 310
 311         if( !reject )
 312         {
 313             t += dt;
 314         }
 315
 316         if( !m_last_step_rejected || boost::numeric::odeint::detail::less_with_sign(new_h, dt, dt) )
 317         {
 318             // limit step size
 319             if( m_max_dt != static_cast<time_type>(0) )
 320             {
 321                 new_h = detail::min_abs(m_max_dt, new_h);
 322             }
 323             m_dt_last = new_h;
 324             dt = new_h;
 325         }
 326
 327         m_last_step_rejected = reject;
 328         m_first = false;
 329
 330         if( reject )
 331             return fail;
 332         else
 333             return success;
 334     }
 335
 336     /** \brief Resets the internal state of the stepper */
 337     void reset()
 338     {
 339         m_first = true;
 340         m_last_step_rejected = false;
 341         // crude estimate of optimal order
 342         m_current_k_opt = 4;
 343         /* no calculation because log10 might not exist for value_type!
 344         const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 );
 345         m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( m_k_max-1 ) , logfact ));
 346         */
 347     }
 348
 349
 350     /* Resizer methods */
 351
 352     template< class StateIn >
 353     void adjust_size( const StateIn &x )
 354     {
 355         resize_m_dxdt( x );
 356         resize_m_xnew( x );
 357         resize_impl( x );
 358         m_midpoint.adjust_size( x );
 359     }
 360
 361
 362 private:
 363
 364     template< class StateIn >
 365     bool resize_m_dxdt( const StateIn &x )
 366     {
 367         return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
 368     }
 369
 370     template< class StateIn >
 371     bool resize_m_xnew( const StateIn &x )
 372     {
 373         return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
 374     }
 375
 376     template< class StateIn >
 377     bool resize_impl( const StateIn &x )
 378     {
 379         bool resized( false );
 380         for( size_t i = 0 ; i < m_k_max ; ++i )
 381             resized |= adjust_size_by_resizeability( m_table[i] , x , typename is_resizeable<state_type>::type() );
 382         resized |= adjust_size_by_resizeability( m_err , x , typename is_resizeable<state_type>::type() );
 383         return resized;
 384     }
 385
 386
 387     template< class System , class StateInOut >
 388     controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
 389     {
 390         typename odeint::unwrap_reference< System >::type &sys = system;
 391         m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_error_bs_type::template resize_m_dxdt< StateInOut > , detail::ref( *this ) , detail::_1 ) );
 392         sys( x , m_dxdt.m_v ,t );
 393         return try_step( system , x , m_dxdt.m_v , t , dt );
 394     }
 395
 396
 397     template< class StateInOut >
 398     void extrapolate( size_t k , state_table_type &table , const value_matrix &coeff , StateInOut &xest )
 399     /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
 400        uses the obtained intermediate results to extrapolate to dt->0
 401     */
 402     {
 403         static const value_type val1 = static_cast< value_type >( 1.0 );
 404         for( int j=k-1 ; j>0 ; --j )
 405         {
 406             m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
 407                                  typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][j] , -coeff[k][j] ) );
 408         }
 409         m_algebra.for_each3( xest , table[0].m_v , xest ,
 410                              typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][0] , -coeff[k][0]) );
 411     }
 412
 413     time_type calc_h_opt( time_type h , value_type error , size_t k ) const
 414     /* calculates the optimal step size for a given error and stage number */
 415     {
 416         BOOST_USING_STD_MIN();
 417         BOOST_USING_STD_MAX();
 418         using std::pow;
 419         value_type expo( 1.0/(2*k+1) );
 420         value_type facmin = m_facmin_table[k];
 421         value_type fac;
 422         if (error == 0.0)
 423             fac=1.0/facmin;
 424         else
 425         {
 426             fac = STEPFAC2 / pow BOOST_PREVENT_MACRO_SUBSTITUTION( error / STEPFAC1 , expo );
 427             fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>(facmin/STEPFAC4) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>(1.0/facmin) , fac ) );
 428         }
 429         return h*fac;
 430     }
 431
 432     controlled_step_result set_k_opt( size_t k , const inv_time_vector &work , const time_vector &h_opt , time_type &dt )
 433     /* calculates the optimal stage number */
 434     {
 435         if( k == 1 )
 436         {
 437             m_current_k_opt = 2;
 438             return success;
 439         }
 440         if( (work[k-1] < KFAC1*work[k]) || (k == m_k_max) )
 441         {   // order decrease
 442             m_current_k_opt = k-1;
 443             dt = h_opt[ m_current_k_opt ];
 444             return success;
 445         }
 446         else if( (work[k] < KFAC2*work[k-1]) || m_last_step_rejected || (k == m_k_max-1) )
 447         {   // same order - also do this if last step got rejected
 448             m_current_k_opt = k;
 449             dt = h_opt[ m_current_k_opt ];
 450             return success;
 451         }
 452         else
 453         {   // order increase - only if last step was not rejected
 454             m_current_k_opt = k+1;
 455             dt = h_opt[ m_current_k_opt-1 ] * m_cost[ m_current_k_opt ] / m_cost[ m_current_k_opt-1 ] ;
 456             return success;
 457         }
 458     }
 459
 460     bool in_convergence_window( size_t k ) const
 461     {
 462         if( (k == m_current_k_opt-1) && !m_last_step_rejected )
 463             return true; // decrease stepsize only if last step was not rejected
 464         return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) );
 465     }
 466
 467     bool should_reject( value_type error , size_t k ) const
 468     {
 469         if( k == m_current_k_opt-1 )
 470         {
 471             const value_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] /
 472                 (m_interval_sequence[0]*m_interval_sequence[0]);
 473             //step will fail, criterion 17.3.17 in NR
 474             return ( error > d*d );
 475         }
 476         else if( k == m_current_k_opt )
 477         {
 478             const value_type d = m_interval_sequence[m_current_k_opt] / m_interval_sequence[0];
 479             return ( error > d*d );
 480         } else
 481             return error > 1.0;
 482     }
 483
 484     default_error_checker< value_type, algebra_type , operations_type > m_error_checker;
 485     modified_midpoint< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint;
 486
 487     bool m_last_step_rejected;
 488     bool m_first;
 489
 490     time_type m_dt_last;
 491     time_type m_t_last;
 492     time_type m_max_dt;
 493
 494     size_t m_current_k_opt;
 495
 496     algebra_type m_algebra;
 497
 498     resizer_type m_dxdt_resizer;
 499     resizer_type m_xnew_resizer;
 500     resizer_type m_resizer;
 501
 502     wrapped_state_type m_xnew;
 503     wrapped_state_type m_err;
 504     wrapped_deriv_type m_dxdt;
 505
 506     int_vector m_interval_sequence; // stores the successive interval counts
 507     value_matrix m_coeff;
 508     int_vector m_cost; // costs for interval count
 509     value_vector m_facmin_table; // for precomputed facmin to save pow calls
 510
 511     state_table_type m_table; // sequence of states for extrapolation
 512
 513     value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
 514 };
 515
 516
 517 /******** DOXYGEN ********/
 518 /**
 519  * \class bulirsch_stoer
 520  * \brief The Bulirsch-Stoer algorithm.
 521  *
 522  * The Bulirsch-Stoer is a controlled stepper that adjusts both step size
 523  * and order of the method. The algorithm uses the modified midpoint and
 524  * a polynomial extrapolation compute the solution.
 525  *
 526  * \tparam State The state type.
 527  * \tparam Value The value type.
 528  * \tparam Deriv The type representing the time derivative of the state.
 529  * \tparam Time The time representing the independent variable - the time.
 530  * \tparam Algebra The algebra type.
 531  * \tparam Operations The operations type.
 532  * \tparam Resizer The resizer policy type.
 533  */
 534
 535     /**
 536      * \fn bulirsch_stoer::bulirsch_stoer( value_type eps_abs , value_type eps_rel , value_type factor_x , value_type factor_dxdt )
 537      * \brief Constructs the bulirsch_stoer class, including initialization of
 538      * the error bounds.
 539      *
 540      * \param eps_abs Absolute tolerance level.
 541      * \param eps_rel Relative tolerance level.
 542      * \param factor_x Factor for the weight of the state.
 543      * \param factor_dxdt Factor for the weight of the derivative.
 544      */
 545
 546     /**
 547      * \fn bulirsch_stoer::try_step( System system , StateInOut &x , time_type &t , time_type &dt )
 548      * \brief Tries to perform one step.
 549      *
 550      * This method tries to do one step with step size dt. If the error estimate
 551      * is to large, the step is rejected and the method returns fail and the
 552      * step size dt is reduced. If the error estimate is acceptably small, the
 553      * step is performed, success is returned and dt might be increased to make
 554      * the steps as large as possible. This method also updates t if a step is
 555      * performed. Also, the internal order of the stepper is adjusted if required.
 556      *
 557      * \param system The system function to solve, hence the r.h.s. of the ODE.
 558      * It must fulfill the Simple System concept.
 559      * \param x The state of the ODE which should be solved. Overwritten if
 560      * the step is successful.
 561      * \param t The value of the time. Updated if the step is successful.
 562      * \param dt The step size. Updated.
 563      * \return success if the step was accepted, fail otherwise.
 564      */
 565
 566     /**
 567      * \fn bulirsch_stoer::try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt )
 568      * \brief Tries to perform one step.
 569      *
 570      * This method tries to do one step with step size dt. If the error estimate
 571      * is to large, the step is rejected and the method returns fail and the
 572      * step size dt is reduced. If the error estimate is acceptably small, the
 573      * step is performed, success is returned and dt might be increased to make
 574      * the steps as large as possible. This method also updates t if a step is
 575      * performed. Also, the internal order of the stepper is adjusted if required.
 576      *
 577      * \param system The system function to solve, hence the r.h.s. of the ODE.
 578      * It must fulfill the Simple System concept.
 579      * \param x The state of the ODE which should be solved. Overwritten if
 580      * the step is successful.
 581      * \param dxdt The derivative of state.
 582      * \param t The value of the time. Updated if the step is successful.
 583      * \param dt The step size. Updated.
 584      * \return success if the step was accepted, fail otherwise.
 585      */
 586
 587     /**
 588      * \fn bulirsch_stoer::try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
 589      * \brief Tries to perform one step.
 590      *
 591      * \note This method is disabled if state_type=time_type to avoid ambiguity.
 592      *
 593      * This method tries to do one step with step size dt. If the error estimate
 594      * is to large, the step is rejected and the method returns fail and the
 595      * step size dt is reduced. If the error estimate is acceptably small, the
 596      * step is performed, success is returned and dt might be increased to make
 597      * the steps as large as possible. This method also updates t if a step is
 598      * performed. Also, the internal order of the stepper is adjusted if required.
 599      *
 600      * \param system The system function to solve, hence the r.h.s. of the ODE.
 601      * It must fulfill the Simple System concept.
 602      * \param in The state of the ODE which should be solved.
 603      * \param t The value of the time. Updated if the step is successful.
 604      * \param out Used to store the result of the step.
 605      * \param dt The step size. Updated.
 606      * \return success if the step was accepted, fail otherwise.
 607      */
 608
 609
 610     /**
 611      * \fn bulirsch_stoer::try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
 612      * \brief Tries to perform one step.
 613      *
 614      * This method tries to do one step with step size dt. If the error estimate
 615      * is to large, the step is rejected and the method returns fail and the
 616      * step size dt is reduced. If the error estimate is acceptably small, the
 617      * step is performed, success is returned and dt might be increased to make
 618      * the steps as large as possible. This method also updates t if a step is
 619      * performed. Also, the internal order of the stepper is adjusted if required.
 620      *
 621      * \param system The system function to solve, hence the r.h.s. of the ODE.
 622      * It must fulfill the Simple System concept.
 623      * \param in The state of the ODE which should be solved.
 624      * \param dxdt The derivative of state.
 625      * \param t The value of the time. Updated if the step is successful.
 626      * \param out Used to store the result of the step.
 627      * \param dt The step size. Updated.
 628      * \return success if the step was accepted, fail otherwise.
 629      */
 630
 631
 632     /**
 633      * \fn bulirsch_stoer::adjust_size( const StateIn &x )
 634      * \brief Adjust the size of all temporaries in the stepper manually.
 635      * \param x A state from which the size of the temporaries to be resized is deduced.
 636      */
 637
 638 }
 639 }
 640 }
 641
 642 #endif // BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED