From 4fb6cf200301f53dc623c0ddd029bd05b4059bf3 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Thu, 4 Aug 2011 09:34:10 +0200 Subject: [PATCH] keine ahnung --- .../numeric/odeint/stepper/bulirsch_stoer.hpp | 406 +++++++++++++++++- 1 file changed, 405 insertions(+), 1 deletion(-) diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp index f840b382..0b64167a 100644 --- a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp +++ b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp @@ -3,7 +3,10 @@ boost/numeric/odeint/stepper/bulirsch_stoer.hpp [begin_description] - Implementaiton of the Burlish-Stoer method. + Implementaiton of the Burlish-Stoer method. As described in + Ernst Hairer, Syvert Paul Nørsett, Gerhard Wanner + Solving Ordinary Differential Equations I. Nonstiff Problems. + Springer Series in Comput. Mathematics, Vol. 8, Springer-Verlag 1987, Second revised edition 1993. [end_description] Copyright 2009-2011 Karsten Ahnert @@ -431,6 +434,407 @@ private: const time_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2; }; + + + + + + + + + + + + +template< +class State , +class Value = double , +class Deriv = State , +class Time = Value , +class Algebra = range_algebra , +class Operations = default_operations , +class Resizer = initially_resizer +> +class bulirsch_stoer_dense_output { + +public: + + typedef State state_type; + typedef Value value_type; + typedef Deriv deriv_type; + typedef Time time_type; + typedef Algebra algebra_type; + typedef Operations operations_type; + typedef Resizer resizer_type; + typedef state_wrapper< state_type > wrapped_state_type; + typedef state_wrapper< deriv_type > wrapped_deriv_type; + typedef controlled_stepper_tag stepper_category; + + typedef bulirsch_stoer< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type; + + typedef std::vector< time_type > value_vector; + typedef std::vector< value_vector > value_matrix; + typedef std::vector< size_t > int_vector; + + bulirsch_stoer_dense_output( + time_type eps_abs = 1E-6 , time_type eps_rel = 1E-6 , + time_type factor_x = 1.0 , time_type factor_dxdt = 1.0 ) + : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ), + m_k_max(8) , + m_safety1(0.25) , m_safety2(0.7), + m_max_dt_factor( 0.1 ) , m_min_step_scale(5E-5) , m_max_step_scale(0.7), + m_last_step_rejected( false ) , m_first( true ) , + m_dt_last( 1.0E30 ) , + m_current_eps( -1.0 ) , + m_error( m_k_max ) , + m_a( m_k_max+1 ) , + m_alpha( m_k_max , value_vector( m_k_max ) ) , + m_interval_sequence( m_k_max+1 ) , + m_coeff( m_k_max+1 ) , + m_cost( m_k_max+1 ) , + m_times( m_k_max ) , + m_table( m_k_max ) , + STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 ) + { + for( unsigned short i = 0; i < m_k_max+1; i++ ) + { + m_interval_sequence[i] = 2 + i*4; + if( i == 0 ) + m_cost[i] = m_interval_sequence[i]; + else + m_cost[i] = m_cost[i-1] + m_interval_sequence[i]; + m_coeff[i].resize(i); + for( size_t k = 0 ; k < i ; ++k ) + { + const time_type r = static_cast< time_type >( m_interval_sequence[i] ) / static_cast< time_type >( m_interval_sequence[k] ); + m_coeff[i][k] = 1.0 / ( r*r - static_cast< time_type >( 1.0 ) ); // coefficients for extrapolation + //std::cout << i << "," << k << " " << m_coeff[i][k] << '\t' ; + } + //std ::cout << std::endl; + // crude estimate of optimal order + const time_type logfact( -log10( std::max( eps_rel , 1.0E-12 ) ) * 0.6 + 0.5 ); + m_current_k_opt = std::max( 1 , std::min( static_cast( m_k_max-1 ) , static_cast( logfact ) )); + //m_current_k_opt = m_k_max - 1; + //std::cout << m_cost[i] << std::endl; + } + + } + + /* + * Version 1 : try_step( sys , x , t , dt ) + * + * The overloads are needed to solve the forwarding problem + */ + template< class System , class StateInOut > + controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt ) + { + return try_step_v1( system , x , t, dt ); + } + + template< class System , class StateInOut > + controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt ) + { + return try_step_v1( system , x , t, dt ); + } + + /* + * Version 2 : try_step( sys , x , dxdt , t , dt ) + * + * this version does not solve the forwarding problem, boost.range can not be used + */ + template< class System , class StateInOut , class DerivIn > + controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt ) + { + m_xnew_resizer.adjust_size( x , boost::bind( &controlled_error_bs_type::template resize_m_xnew< StateInOut > , boost::ref( *this ) , _1 ) ); + controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt ); + if( ( res == success_step_size_increased ) || ( res == success_step_size_unchanged ) ) + { + boost::numeric::odeint::copy( m_xnew.m_v , x ); + } + return res; + } + + /* + * Version 3 : try_step( sys , in , t , out , dt ) + * + * this version does not solve the forwarding problem, boost.range can not be used + */ + template< class System , class StateIn , class StateOut > + controlled_step_result try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt ) + { + typename boost::unwrap_reference< System >::type &sys = system; + m_dxdt_resizer.adjust_size( in , boost::bind( &controlled_error_bs_type::template resize_m_dxdt< StateIn > , boost::ref( *this ) , _1 ) ); + sys( in , m_dxdt.m_v , t ); + return try_step( system , in , m_dxdt.m_v , t , out , dt ); + } + + + template< class System , class StateIn , class DerivIn , class StateOut > + controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt ) + { + static const time_type val1( static_cast< time_type >( 1.0 ) ); + + typename boost::unwrap_reference< System >::type &sys = system; + if( m_resizer.adjust_size( in , boost::bind( &controlled_error_bs_type::template resize< StateIn > , boost::ref( *this ) , _1 ) ) ) + reset(); // system resized -> reset + if( dt != m_dt_last ) + reset(); // step size changed from outside -> reset + + bool reject( true ); + m_dt_last = dt; + + value_vector h_opt( m_k_max+1 ); + value_vector work( m_k_max+1 ); + + size_t k_final = 0; + + //std::cout << "t=" << t <<", dt=" << dt << ", k_opt=" << m_current_k_opt << std::endl; + + for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ ) + { + //std::cout << "k=" << k <<": " << std::endl; + m_midpoint.set_steps( m_interval_sequence[k] ); + if( k == 0 ) + { + m_midpoint.do_step( sys , in , dxdt , t , out , dt ); + } + else + { + m_midpoint.do_step( sys , in , dxdt , t , m_table[k-1].m_v , dt ); + extrapolate( k , out ); + // get error estimate + m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v , + typename operations_type::template scale_sum2< time_type , time_type >( val1 , -val1 ) ); + const time_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt ); + h_opt[k] = calc_h_opt( dt , error , k ); + work[k] = m_cost[k]/h_opt[k]; + //std::cout << '\t' << "h_opt=" << h_opt[k] << ", work=" << work[k] << std::endl; + //std::cout << '\t' << "error: " << error << std::endl; + if( m_first && (error <= static_cast< time_type >( 1.0 )) ) + { // this is the first step, convergence does not have to be in order window + //std::cout << '\t' << "convergence in first step" << std::endl; + reject = false; + k_final = k; + break; // leave k-loop + } + if( in_convergence_window( k ) ) + { + k_final = k; + if( error <= 1.0 ) // convergence in convergence window + { + reject = false; + break; // we are good, leave k-loop + } + else + { + if( should_reject( error , k ) ) + { + reject = true; + break; + } + } + } + } + } + + //calculate optimal order and stepsize + controlled_step_result result = set_k_opt( k_final , work , h_opt , dt ); + + m_last_step_rejected = reject; + if( reject ) + result = step_size_decreased; + else + t += m_dt_last; + + m_dt_last = dt; + + m_first = false; + + return result; + } + + void reset() + { + m_first = true; + m_last_step_rejected = false; + } + + + /* Resizer methods */ + + template< class StateIn > + bool resize_m_dxdt( const StateIn &x ) + { + return adjust_size_by_resizeability( m_dxdt , x , typename wrapped_deriv_type::is_resizeable() ); + } + + template< class StateIn > + bool resize_m_xnew( const StateIn &x ) + { + return adjust_size_by_resizeability( m_xnew , x , typename wrapped_state_type::is_resizeable() ); + } + + template< class StateIn > + bool resize( const StateIn &x ) + { + bool resized( false ); + for( size_t i = 0 ; i < m_k_max ; ++i ) + resized |= adjust_size_by_resizeability( m_table[i] , x , typename wrapped_state_type::is_resizeable() ); + resized |= adjust_size_by_resizeability( m_err , x , typename wrapped_state_type::is_resizeable() ); + return resized; + } + + template< class StateIn > + void adjust_size( const StateIn &x ) + { + resize_m_dxdt( x ); + resize_m_xnew( x ); + resize( x ); + m_midpoint.adjust_size(); + } + + +private: + + template< class System , class StateInOut > + controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt ) + { + typename boost::unwrap_reference< System >::type &sys = system; + m_dxdt_resizer.adjust_size( x , boost::bind( &controlled_error_bs_type::template resize_m_dxdt< StateInOut > , boost::ref( *this ) , _1 ) ); + sys( x , m_dxdt.m_v ,t ); + return try_step( system , x , m_dxdt.m_v , t , dt ); + } + + + template< class StateInOut > + void extrapolate( const size_t k , StateInOut &xest ) + //polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf + { + //std::cout << "extrapolate k=" << k << ":" << std::endl; + static const time_type val1 = static_cast< time_type >( 1.0 ); + for( int j=k-1 ; j>0 ; --j ) + { + //std::cout << '\t' << m_coeff[k][j]; + m_algebra.for_each3( m_table[j-1].m_v , m_table[j].m_v , m_table[j-1].m_v , + typename operations_type::template scale_sum2< time_type , time_type >( val1 + m_coeff[k][j] , -m_coeff[k][j] ) ); + } + //std::cout << std::endl << m_coeff[k][0] << std::endl; + m_algebra.for_each3( xest , m_table[0].m_v , xest , + typename operations_type::template scale_sum2< time_type , time_type >( val1 + m_coeff[k][0] , -m_coeff[k][0]) ); + } + + time_type calc_h_opt( const time_type h , const value_type error , const size_t k ) const + { + time_type expo=1.0/(2*k+1); + time_type facmin = std::pow( STEPFAC3 , expo ); + time_type fac; + if (error == 0.0) + fac=1.0/facmin; + else + { + fac = STEPFAC2 / std::pow( error / STEPFAC1 , expo ); + fac = std::max( facmin/STEPFAC4 , std::min( 1.0/facmin , fac ) ); + } + return std::abs(h*fac); + } + + controlled_step_result set_k_opt( const size_t k , const value_vector &work , const value_vector &h_opt , time_type &dt ) + { + //std::cout << "finding k_opt..." << std::endl; + if( k == 1 ) + { + m_current_k_opt = 2; + //dt = h_opt[ m_current_k_opt-1 ] * m_cost[ m_current_k_opt ] / m_cost[ m_current_k_opt-1 ] ; + return success_step_size_increased; + } + if( (work[k-1] < KFAC1*work[k]) || (k == m_k_max) ) + { // order decrease + m_current_k_opt = k-1; + dt = h_opt[ m_current_k_opt ]; + return success_step_size_increased; + } + else if( (work[k] < KFAC2*work[k-1]) || m_last_step_rejected || (k == m_k_max-1) ) + { // same order - also do this if last step got rejected + m_current_k_opt = k; + dt = h_opt[ m_current_k_opt ]; + return success_step_size_unchanged; + } + else + { // order increase - only if last step was not rejected + m_current_k_opt = k+1; + dt = h_opt[ m_current_k_opt-1 ] * m_cost[ m_current_k_opt ] / m_cost[ m_current_k_opt-1 ] ; + return success_step_size_increased; + } + } + + bool in_convergence_window( const size_t k ) const + { + if( (k == m_current_k_opt-1) && !m_last_step_rejected ) + return true; // decrease stepsize only if last step was not rejected + return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) ); + } + + bool should_reject( const time_type error , const size_t k ) const + { + if( (k == m_current_k_opt-1) ) + { + const time_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] / + (m_interval_sequence[0]*m_interval_sequence[0]); + //step will fail, criterion 17.3.17 in NR + return ( error > d*d ); + } + else if( k == m_current_k_opt ) + { + const time_type d = m_interval_sequence[m_current_k_opt] / m_interval_sequence[0]; + return ( error > d*d ); + } else + return error > 1.0; + } + + default_error_checker< value_type, algebra_type , operations_type > m_error_checker; + modified_midpoint< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint; + + const size_t m_k_max; + + const time_type m_safety1; + const time_type m_safety2; + const time_type m_max_dt_factor; + const time_type m_min_step_scale; + const time_type m_max_step_scale; + + bool m_last_step_rejected; + bool m_first; + + time_type m_dt_last; + time_type m_t_last; + time_type m_current_eps; + + size_t m_current_k_max; + size_t m_current_k_opt; + + algebra_type m_algebra; + + resizer_type m_dxdt_resizer; + resizer_type m_xnew_resizer; + resizer_type m_resizer; + + wrapped_state_type m_xnew; + wrapped_state_type m_err; + wrapped_deriv_type m_dxdt; + + value_vector m_error; // errors of repeated midpoint steps and extrapolations + value_vector m_a; // stores the work (number of f calls) required for the orders + value_matrix m_alpha; // stores convergence factor for stepsize adjustment + int_vector m_interval_sequence; // stores the successive interval counts + value_matrix m_coeff; + int_vector m_cost; // costs for interval count + + value_vector m_times; + std::vector< wrapped_state_type > m_table; // sequence of states for extrapolation + + const time_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2; +}; + } } }