Chris@16: /*
Chris@16:   [auto_generated]
Chris@16:   boost/numeric/odeint/stepper/bulirsch_stoer.hpp
Chris@16: 
Chris@16:   [begin_description]
Chris@16:   Implementation of the Burlish-Stoer method. As described in
Chris@16:   Ernst Hairer, Syvert Paul Norsett, Gerhard Wanner
Chris@16:   Solving Ordinary Differential Equations I. Nonstiff Problems.
Chris@16:   Springer Series in Comput. Mathematics, Vol. 8, Springer-Verlag 1987, Second revised edition 1993.
Chris@16:   [end_description]
Chris@16: 
Chris@101:   Copyright 2011-2013 Mario Mulansky
Chris@101:   Copyright 2011-2013 Karsten Ahnert
Chris@101:   Copyright 2012 Christoph Koke
Chris@16: 
Chris@16:   Distributed under the Boost Software License, Version 1.0.
Chris@16:   (See accompanying file LICENSE_1_0.txt or
Chris@16:   copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: */
Chris@16: 
Chris@16: 
Chris@16: #ifndef BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED
Chris@16: #define BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED
Chris@16: 
Chris@16: 
Chris@16: #include <iostream>
Chris@16: 
Chris@16: #include <algorithm>
Chris@16: 
Chris@16: #include <boost/config.hpp> // for min/max guidelines
Chris@16: 
Chris@16: #include <boost/numeric/odeint/util/bind.hpp>
Chris@16: #include <boost/numeric/odeint/util/unwrap_reference.hpp>
Chris@16: 
Chris@16: #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
Chris@16: #include <boost/numeric/odeint/stepper/modified_midpoint.hpp>
Chris@16: #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
Chris@16: #include <boost/numeric/odeint/algebra/range_algebra.hpp>
Chris@16: #include <boost/numeric/odeint/algebra/default_operations.hpp>
Chris@101: #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
Chris@101: #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
Chris@16: 
Chris@16: #include <boost/numeric/odeint/util/state_wrapper.hpp>
Chris@16: #include <boost/numeric/odeint/util/is_resizeable.hpp>
Chris@16: #include <boost/numeric/odeint/util/resizer.hpp>
Chris@16: #include <boost/numeric/odeint/util/unit_helper.hpp>
Chris@16: #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
Chris@16: 
Chris@16: namespace boost {
Chris@16: namespace numeric {
Chris@16: namespace odeint {
Chris@16: 
Chris@16: template<
Chris@16:     class State ,
Chris@16:     class Value = double ,
Chris@16:     class Deriv = State ,
Chris@16:     class Time = Value ,
Chris@101:     class Algebra = typename algebra_dispatcher< State >::algebra_type ,
Chris@101:     class Operations = typename operations_dispatcher< State >::operations_type ,
Chris@16:     class Resizer = initially_resizer
Chris@16:     >
Chris@16: class bulirsch_stoer {
Chris@16: 
Chris@16: public:
Chris@16: 
Chris@16:     typedef State state_type;
Chris@16:     typedef Value value_type;
Chris@16:     typedef Deriv deriv_type;
Chris@16:     typedef Time time_type;
Chris@16:     typedef Algebra algebra_type;
Chris@16:     typedef Operations operations_type;
Chris@16:     typedef Resizer resizer_type;
Chris@16: #ifndef DOXYGEN_SKIP
Chris@16:     typedef state_wrapper< state_type > wrapped_state_type;
Chris@16:     typedef state_wrapper< deriv_type > wrapped_deriv_type;
Chris@16:     typedef controlled_stepper_tag stepper_category;
Chris@16: 
Chris@16:     typedef bulirsch_stoer< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type;
Chris@16: 
Chris@16:     typedef typename inverse_time< time_type >::type inv_time_type;
Chris@16: 
Chris@16:     typedef std::vector< value_type > value_vector;
Chris@16:     typedef std::vector< time_type > time_vector;
Chris@16:     typedef std::vector< inv_time_type > inv_time_vector;  //should be 1/time_type for boost.units
Chris@16:     typedef std::vector< value_vector > value_matrix;
Chris@16:     typedef std::vector< size_t > int_vector;
Chris@16:     typedef std::vector< wrapped_state_type > state_table_type;
Chris@16: #endif //DOXYGEN_SKIP
Chris@16:     const static size_t m_k_max = 8;
Chris@16: 
Chris@16:     bulirsch_stoer(
Chris@16:         value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
Chris@16:         value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 )
Chris@16:         : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , m_midpoint() ,
Chris@16:           m_last_step_rejected( false ) , m_first( true ) ,
Chris@16:           m_interval_sequence( m_k_max+1 ) ,
Chris@16:           m_coeff( m_k_max+1 ) ,
Chris@16:           m_cost( m_k_max+1 ) ,
Chris@16:           m_table( m_k_max ) ,
Chris@16:           STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 )
Chris@16:     {
Chris@16:         BOOST_USING_STD_MIN();
Chris@16:         BOOST_USING_STD_MAX();
Chris@16:         /* initialize sequence of stage numbers and work */
Chris@16:         for( unsigned short i = 0; i < m_k_max+1; i++ )
Chris@16:         {
Chris@16:             m_interval_sequence[i] = 2 * (i+1);
Chris@16:             if( i == 0 )
Chris@16:                 m_cost[i] = m_interval_sequence[i];
Chris@16:             else
Chris@16:                 m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
Chris@16:             m_coeff[i].resize(i);
Chris@16:             for( size_t k = 0 ; k < i ; ++k  )
Chris@16:             {
Chris@16:                 const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
Chris@16:                 m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
Chris@16:             }
Chris@16: 
Chris@16:             // crude estimate of optimal order
Chris@16: 
Chris@16:             m_current_k_opt = 4;
Chris@16:             /* no calculation because log10 might not exist for value_type!
Chris@16:             const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 );
Chris@16:             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 ));
Chris@16:             */
Chris@16:         }
Chris@16: 
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16:     /*
Chris@16:      * Version 1 : try_step( sys , x , t , dt )
Chris@16:      *
Chris@16:      * The overloads are needed to solve the forwarding problem
Chris@16:      */
Chris@16:     template< class System , class StateInOut >
Chris@16:     controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
Chris@16:     {
Chris@16:         return try_step_v1( system , x , t, dt );
Chris@16:     }
Chris@16: 
Chris@16:     /**
Chris@16:      * \brief Second version to solve the forwarding problem, can be used with Boost.Range as StateInOut.
Chris@16:      */
Chris@16:     template< class System , class StateInOut >
Chris@16:     controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
Chris@16:     {
Chris@16:         return try_step_v1( system , x , t, dt );
Chris@16:     }
Chris@16: 
Chris@16:     /*
Chris@16:      * Version 2 : try_step( sys , x , dxdt , t , dt )
Chris@16:      *
Chris@16:      * this version does not solve the forwarding problem, boost.range can not be used
Chris@16:      */
Chris@16:     template< class System , class StateInOut , class DerivIn >
Chris@16:     controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt )
Chris@16:     {
Chris@16:         m_xnew_resizer.adjust_size( x , detail::bind( &controlled_error_bs_type::template resize_m_xnew< StateInOut > , detail::ref( *this ) , detail::_1 ) );
Chris@16:         controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt );
Chris@16:         if( res == success )
Chris@16:         {
Chris@16:             boost::numeric::odeint::copy( m_xnew.m_v , x );
Chris@16:         }
Chris@16:         return res;
Chris@16:     }
Chris@16: 
Chris@16:     /*
Chris@16:      * Version 3 : try_step( sys , in , t , out , dt )
Chris@16:      *
Chris@16:      * this version does not solve the forwarding problem, boost.range can not be used
Chris@16:      */
Chris@16:     template< class System , class StateIn , class StateOut >
Chris@16:     typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type
Chris@16:     try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
Chris@16:     {
Chris@16:         typename odeint::unwrap_reference< System >::type &sys = system;
Chris@16:         m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_m_dxdt< StateIn > , detail::ref( *this ) , detail::_1 ) );
Chris@16:         sys( in , m_dxdt.m_v , t );
Chris@16:         return try_step( system , in , m_dxdt.m_v , t , out , dt );
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16:     /*
Chris@16:      * Full version : try_step( sys , in , dxdt_in , t , out , dt )
Chris@16:      *
Chris@16:      * contains the actual implementation
Chris@16:      */
Chris@16:     template< class System , class StateIn , class DerivIn , class StateOut >
Chris@16:     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
Chris@16:     {
Chris@16:         BOOST_USING_STD_MIN();
Chris@16:         BOOST_USING_STD_MAX();
Chris@16: 
Chris@16:         static const value_type val1( 1.0 );
Chris@16: 
Chris@16:         if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) )
Chris@16:         {
Chris@16:             reset(); // system resized -> reset
Chris@16:         }
Chris@16: 
Chris@16:         if( dt != m_dt_last )
Chris@16:         {
Chris@16:             reset(); // step size changed from outside -> reset
Chris@16:         }
Chris@16: 
Chris@16:         bool reject( true );
Chris@16: 
Chris@16:         time_vector h_opt( m_k_max+1 );
Chris@16:         inv_time_vector work( m_k_max+1 );
Chris@16: 
Chris@16:         time_type new_h = dt;
Chris@16: 
Chris@16:         /* m_current_k_opt is the estimated current optimal stage number */
Chris@16:         for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ )
Chris@16:         {
Chris@16:             /* the stage counts are stored in m_interval_sequence */
Chris@16:             m_midpoint.set_steps( m_interval_sequence[k] );
Chris@16:             if( k == 0 )
Chris@16:             {
Chris@101:                 m_midpoint.do_step( system , in , dxdt , t , out , dt );
Chris@16:                 /* the first step, nothing more to do */
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@101:                 m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt );
Chris@16:                 extrapolate( k , m_table , m_coeff , out );
Chris@16:                 // get error estimate
Chris@16:                 m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v ,
Chris@16:                                      typename operations_type::template scale_sum2< value_type , value_type >( val1 , -val1 ) );
Chris@16:                 const value_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt );
Chris@16:                 h_opt[k] = calc_h_opt( dt , error , k );
Chris@16:                 work[k] = static_cast<value_type>( m_cost[k] ) / h_opt[k];
Chris@16: 
Chris@16:                 if( (k == m_current_k_opt-1) || m_first )
Chris@16:                 { // convergence before k_opt ?
Chris@16:                     if( error < 1.0 )
Chris@16:                     {
Chris@16:                         //convergence
Chris@16:                         reject = false;
Chris@16:                         if( (work[k] < KFAC2*work[k-1]) || (m_current_k_opt <= 2) )
Chris@16:                         {
Chris@16:                             // leave order as is (except we were in first round)
Chris@16:                             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 ) );
Chris@16:                             new_h = h_opt[k];
Chris@16:                             new_h *= static_cast<value_type>( m_cost[k+1] ) / static_cast<value_type>( m_cost[k] );
Chris@16:                         } else {
Chris@16:                             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) ) );
Chris@16:                             new_h = h_opt[k];
Chris@16:                         }
Chris@16:                         break;
Chris@16:                     }
Chris@16:                     else if( should_reject( error , k ) && !m_first )
Chris@16:                     {
Chris@16:                         reject = true;
Chris@16:                         new_h = h_opt[k];
Chris@16:                         break;
Chris@16:                     }
Chris@16:                 }
Chris@16:                 if( k == m_current_k_opt )
Chris@16:                 { // convergence at k_opt ?
Chris@16:                     if( error < 1.0 )
Chris@16:                     {
Chris@16:                         //convergence
Chris@16:                         reject = false;
Chris@16:                         if( (work[k-1] < KFAC2*work[k]) )
Chris@16:                         {
Chris@16:                             m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
Chris@16:                             new_h = h_opt[m_current_k_opt];
Chris@16:                         }
Chris@16:                         else if( (work[k] < KFAC2*work[k-1]) && !m_last_step_rejected )
Chris@16:                         {
Chris@16:                             m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max-1) , static_cast<int>(m_current_k_opt)+1 );
Chris@16:                             new_h = h_opt[k];
Chris@16:                             new_h *= m_cost[m_current_k_opt]/m_cost[k];
Chris@16:                         } else
Chris@16:                             new_h = h_opt[m_current_k_opt];
Chris@16:                         break;
Chris@16:                     }
Chris@16:                     else if( should_reject( error , k ) )
Chris@16:                     {
Chris@16:                         reject = true;
Chris@16:                         new_h = h_opt[m_current_k_opt];
Chris@16:                         break;
Chris@16:                     }
Chris@16:                 }
Chris@16:                 if( k == m_current_k_opt+1 )
Chris@16:                 { // convergence at k_opt+1 ?
Chris@16:                     if( error < 1.0 )
Chris@16:                     {   //convergence
Chris@16:                         reject = false;
Chris@16:                         if( work[k-2] < KFAC2*work[k-1] )
Chris@16:                             m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
Chris@16:                         if( (work[k] < KFAC2*work[m_current_k_opt]) && !m_last_step_rejected )
Chris@16:                             m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , static_cast<int>(k) );
Chris@16:                         new_h = h_opt[m_current_k_opt];
Chris@16:                     } else
Chris@16:                     {
Chris@16:                         reject = true;
Chris@16:                         new_h = h_opt[m_current_k_opt];
Chris@16:                     }
Chris@16:                     break;
Chris@16:                 }
Chris@16:             }
Chris@16:         }
Chris@16: 
Chris@16:         if( !reject )
Chris@16:         {
Chris@16:             t += dt;
Chris@16:         }
Chris@16: 
Chris@16:         if( !m_last_step_rejected || boost::numeric::odeint::detail::less_with_sign(new_h, dt, dt) )
Chris@16:         {
Chris@16:             m_dt_last = new_h;
Chris@16:             dt = new_h;
Chris@16:         }
Chris@16: 
Chris@16:         m_last_step_rejected = reject;
Chris@16:         m_first = false;
Chris@16: 
Chris@16:         if( reject )
Chris@16:             return fail;
Chris@16:         else
Chris@16:             return success;
Chris@16:     }
Chris@16: 
Chris@16:     /** \brief Resets the internal state of the stepper */
Chris@16:     void reset()
Chris@16:     {
Chris@16:         m_first = true;
Chris@16:         m_last_step_rejected = false;
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16:     /* Resizer methods */
Chris@16: 
Chris@16:     template< class StateIn >
Chris@16:     void adjust_size( const StateIn &x )
Chris@16:     {
Chris@16:         resize_m_dxdt( x );
Chris@16:         resize_m_xnew( x );
Chris@16:         resize_impl( x );
Chris@101:         m_midpoint.adjust_size( x );
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16: private:
Chris@16: 
Chris@16:     template< class StateIn >
Chris@16:     bool resize_m_dxdt( const StateIn &x )
Chris@16:     {
Chris@16:         return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
Chris@16:     }
Chris@16: 
Chris@16:     template< class StateIn >
Chris@16:     bool resize_m_xnew( const StateIn &x )
Chris@16:     {
Chris@16:         return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
Chris@16:     }
Chris@16: 
Chris@16:     template< class StateIn >
Chris@16:     bool resize_impl( const StateIn &x )
Chris@16:     {
Chris@16:         bool resized( false );
Chris@16:         for( size_t i = 0 ; i < m_k_max ; ++i )
Chris@16:             resized |= adjust_size_by_resizeability( m_table[i] , x , typename is_resizeable<state_type>::type() );
Chris@16:         resized |= adjust_size_by_resizeability( m_err , x , typename is_resizeable<state_type>::type() );
Chris@16:         return resized;
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16:     template< class System , class StateInOut >
Chris@16:     controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
Chris@16:     {
Chris@16:         typename odeint::unwrap_reference< System >::type &sys = system;
Chris@16:         m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_error_bs_type::template resize_m_dxdt< StateInOut > , detail::ref( *this ) , detail::_1 ) );
Chris@16:         sys( x , m_dxdt.m_v ,t );
Chris@16:         return try_step( system , x , m_dxdt.m_v , t , dt );
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16:     template< class StateInOut >
Chris@16:     void extrapolate( size_t k , state_table_type &table , const value_matrix &coeff , StateInOut &xest )
Chris@16:     /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
Chris@16:        uses the obtained intermediate results to extrapolate to dt->0 
Chris@16:     */
Chris@16:     {
Chris@16:         static const value_type val1 = static_cast< value_type >( 1.0 );
Chris@16:         for( int j=k-1 ; j>0 ; --j )
Chris@16:         {
Chris@16:             m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
Chris@16:                                  typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][j] , -coeff[k][j] ) );
Chris@16:         }
Chris@16:         m_algebra.for_each3( xest , table[0].m_v , xest ,
Chris@16:                              typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][0] , -coeff[k][0]) );
Chris@16:     }
Chris@16: 
Chris@16:     time_type calc_h_opt( time_type h , value_type error , size_t k ) const
Chris@16:     /* calculates the optimal step size for a given error and stage number */
Chris@16:     {
Chris@16:         BOOST_USING_STD_MIN();
Chris@16:         BOOST_USING_STD_MAX();
Chris@16:         using std::pow;
Chris@16:         value_type expo( 1.0/(2*k+1) );
Chris@16:         value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
Chris@16:         value_type fac;
Chris@16:         if (error == 0.0)
Chris@16:             fac=1.0/facmin;
Chris@16:         else
Chris@16:         {
Chris@16:             fac = STEPFAC2 / pow BOOST_PREVENT_MACRO_SUBSTITUTION( error / STEPFAC1 , expo );
Chris@101:             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 ) );
Chris@16:         }
Chris@16:         return h*fac;
Chris@16:     }
Chris@16: 
Chris@16:     controlled_step_result set_k_opt( size_t k , const inv_time_vector &work , const time_vector &h_opt , time_type &dt )
Chris@16:     /* calculates the optimal stage number */
Chris@16:     {
Chris@16:         if( k == 1 )
Chris@16:         {
Chris@16:             m_current_k_opt = 2;
Chris@16:             return success;
Chris@16:         }
Chris@16:         if( (work[k-1] < KFAC1*work[k]) || (k == m_k_max) )
Chris@16:         {   // order decrease
Chris@16:             m_current_k_opt = k-1;
Chris@16:             dt = h_opt[ m_current_k_opt ];
Chris@16:             return success;
Chris@16:         }
Chris@16:         else if( (work[k] < KFAC2*work[k-1]) || m_last_step_rejected || (k == m_k_max-1) )
Chris@16:         {   // same order - also do this if last step got rejected
Chris@16:             m_current_k_opt = k;
Chris@16:             dt = h_opt[ m_current_k_opt ];
Chris@16:             return success;
Chris@16:         }
Chris@16:         else
Chris@16:         {   // order increase - only if last step was not rejected
Chris@16:             m_current_k_opt = k+1;
Chris@16:             dt = h_opt[ m_current_k_opt-1 ] * m_cost[ m_current_k_opt ] / m_cost[ m_current_k_opt-1 ] ;
Chris@16:             return success;
Chris@16:         }
Chris@16:     }
Chris@16: 
Chris@16:     bool in_convergence_window( size_t k ) const
Chris@16:     {
Chris@16:         if( (k == m_current_k_opt-1) && !m_last_step_rejected )
Chris@16:             return true; // decrease stepsize only if last step was not rejected
Chris@16:         return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) );
Chris@16:     }
Chris@16: 
Chris@16:     bool should_reject( value_type error , size_t k ) const
Chris@16:     {
Chris@16:         if( k == m_current_k_opt-1 )
Chris@16:         {
Chris@16:             const value_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] /
Chris@16:                 (m_interval_sequence[0]*m_interval_sequence[0]);
Chris@16:             //step will fail, criterion 17.3.17 in NR
Chris@16:             return ( error > d*d );
Chris@16:         }
Chris@16:         else if( k == m_current_k_opt )
Chris@16:         {
Chris@16:             const value_type d = m_interval_sequence[m_current_k_opt] / m_interval_sequence[0];
Chris@16:             return ( error > d*d );
Chris@16:         } else
Chris@16:             return error > 1.0;
Chris@16:     }
Chris@16: 
Chris@16:     default_error_checker< value_type, algebra_type , operations_type > m_error_checker;
Chris@16:     modified_midpoint< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint;
Chris@16: 
Chris@16:     bool m_last_step_rejected;
Chris@16:     bool m_first;
Chris@16: 
Chris@16:     time_type m_dt_last;
Chris@16:     time_type m_t_last;
Chris@16: 
Chris@16:     size_t m_current_k_opt;
Chris@16: 
Chris@16:     algebra_type m_algebra;
Chris@16: 
Chris@16:     resizer_type m_dxdt_resizer;
Chris@16:     resizer_type m_xnew_resizer;
Chris@16:     resizer_type m_resizer;
Chris@16: 
Chris@16:     wrapped_state_type m_xnew;
Chris@16:     wrapped_state_type m_err;
Chris@16:     wrapped_deriv_type m_dxdt;
Chris@16: 
Chris@16:     int_vector m_interval_sequence; // stores the successive interval counts
Chris@16:     value_matrix m_coeff;
Chris@16:     int_vector m_cost; // costs for interval count
Chris@16: 
Chris@16:     state_table_type m_table; // sequence of states for extrapolation
Chris@16: 
Chris@16:     const value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: /******** DOXYGEN ********/
Chris@16: /**
Chris@16:  * \class bulirsch_stoer
Chris@16:  * \brief The Bulirsch-Stoer algorithm.
Chris@16:  * 
Chris@16:  * The Bulirsch-Stoer is a controlled stepper that adjusts both step size
Chris@16:  * and order of the method. The algorithm uses the modified midpoint and
Chris@16:  * a polynomial extrapolation compute the solution.
Chris@16:  *
Chris@16:  * \tparam State The state type.
Chris@16:  * \tparam Value The value type.
Chris@16:  * \tparam Deriv The type representing the time derivative of the state.
Chris@16:  * \tparam Time The time representing the independent variable - the time.
Chris@16:  * \tparam Algebra The algebra type.
Chris@16:  * \tparam Operations The operations type.
Chris@16:  * \tparam Resizer The resizer policy type.
Chris@16:  */
Chris@16: 
Chris@16:     /**
Chris@16:      * \fn bulirsch_stoer::bulirsch_stoer( value_type eps_abs , value_type eps_rel , value_type factor_x , value_type factor_dxdt )
Chris@16:      * \brief Constructs the bulirsch_stoer class, including initialization of 
Chris@16:      * the error bounds.
Chris@16:      *
Chris@16:      * \param eps_abs Absolute tolerance level.
Chris@16:      * \param eps_rel Relative tolerance level.
Chris@16:      * \param factor_x Factor for the weight of the state.
Chris@16:      * \param factor_dxdt Factor for the weight of the derivative.
Chris@16:      */
Chris@16: 
Chris@16:     /**
Chris@16:      * \fn bulirsch_stoer::try_step( System system , StateInOut &x , time_type &t , time_type &dt )
Chris@16:      * \brief Tries to perform one step.
Chris@16:      *
Chris@16:      * This method tries to do one step with step size dt. If the error estimate
Chris@16:      * is to large, the step is rejected and the method returns fail and the 
Chris@16:      * step size dt is reduced. If the error estimate is acceptably small, the
Chris@16:      * step is performed, success is returned and dt might be increased to make 
Chris@16:      * the steps as large as possible. This method also updates t if a step is
Chris@16:      * performed. Also, the internal order of the stepper is adjusted if required.
Chris@16:      *
Chris@16:      * \param system The system function to solve, hence the r.h.s. of the ODE. 
Chris@16:      * It must fulfill the Simple System concept.
Chris@16:      * \param x The state of the ODE which should be solved. Overwritten if 
Chris@16:      * the step is successful.
Chris@16:      * \param t The value of the time. Updated if the step is successful.
Chris@16:      * \param dt The step size. Updated.
Chris@16:      * \return success if the step was accepted, fail otherwise.
Chris@16:      */
Chris@16: 
Chris@16:     /**
Chris@16:      * \fn bulirsch_stoer::try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt )
Chris@16:      * \brief Tries to perform one step.
Chris@16:      *
Chris@16:      * This method tries to do one step with step size dt. If the error estimate
Chris@16:      * is to large, the step is rejected and the method returns fail and the 
Chris@16:      * step size dt is reduced. If the error estimate is acceptably small, the
Chris@16:      * step is performed, success is returned and dt might be increased to make 
Chris@16:      * the steps as large as possible. This method also updates t if a step is
Chris@16:      * performed. Also, the internal order of the stepper is adjusted if required.
Chris@16:      *
Chris@16:      * \param system The system function to solve, hence the r.h.s. of the ODE. 
Chris@16:      * It must fulfill the Simple System concept.
Chris@16:      * \param x The state of the ODE which should be solved. Overwritten if 
Chris@16:      * the step is successful.
Chris@16:      * \param dxdt The derivative of state.
Chris@16:      * \param t The value of the time. Updated if the step is successful.
Chris@16:      * \param dt The step size. Updated.
Chris@16:      * \return success if the step was accepted, fail otherwise.
Chris@16:      */
Chris@16: 
Chris@16:     /**
Chris@16:      * \fn bulirsch_stoer::try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
Chris@16:      * \brief Tries to perform one step.
Chris@16:      *
Chris@16:      * \note This method is disabled if state_type=time_type to avoid ambiguity.
Chris@16:      *
Chris@16:      * This method tries to do one step with step size dt. If the error estimate
Chris@16:      * is to large, the step is rejected and the method returns fail and the 
Chris@16:      * step size dt is reduced. If the error estimate is acceptably small, the
Chris@16:      * step is performed, success is returned and dt might be increased to make 
Chris@16:      * the steps as large as possible. This method also updates t if a step is
Chris@16:      * performed. Also, the internal order of the stepper is adjusted if required.
Chris@16:      *
Chris@16:      * \param system The system function to solve, hence the r.h.s. of the ODE. 
Chris@16:      * It must fulfill the Simple System concept.
Chris@16:      * \param in The state of the ODE which should be solved.
Chris@16:      * \param t The value of the time. Updated if the step is successful.
Chris@16:      * \param out Used to store the result of the step.
Chris@16:      * \param dt The step size. Updated.
Chris@16:      * \return success if the step was accepted, fail otherwise.
Chris@16:      */
Chris@16: 
Chris@16: 
Chris@16:     /**
Chris@16:      * \fn bulirsch_stoer::try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
Chris@16:      * \brief Tries to perform one step.
Chris@16:      *
Chris@16:      * This method tries to do one step with step size dt. If the error estimate
Chris@16:      * is to large, the step is rejected and the method returns fail and the 
Chris@16:      * step size dt is reduced. If the error estimate is acceptably small, the
Chris@16:      * step is performed, success is returned and dt might be increased to make 
Chris@16:      * the steps as large as possible. This method also updates t if a step is
Chris@16:      * performed. Also, the internal order of the stepper is adjusted if required.
Chris@16:      *
Chris@16:      * \param system The system function to solve, hence the r.h.s. of the ODE. 
Chris@16:      * It must fulfill the Simple System concept.
Chris@16:      * \param in The state of the ODE which should be solved.
Chris@16:      * \param dxdt The derivative of state.
Chris@16:      * \param t The value of the time. Updated if the step is successful.
Chris@16:      * \param out Used to store the result of the step.
Chris@16:      * \param dt The step size. Updated.
Chris@16:      * \return success if the step was accepted, fail otherwise.
Chris@16:      */
Chris@16: 
Chris@16: 
Chris@16:     /**
Chris@16:      * \fn bulirsch_stoer::adjust_size( const StateIn &x )
Chris@16:      * \brief Adjust the size of all temporaries in the stepper manually.
Chris@16:      * \param x A state from which the size of the temporaries to be resized is deduced.
Chris@16:      */
Chris@16: 
Chris@16: }
Chris@16: }
Chris@16: }
Chris@16: 
Chris@16: #endif // BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED