3 boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
6 Implementaiton of the Burlish-Stoer method with dense output
9 Copyright 2011-2013 Mario Mulansky
10 Copyright 2011-2013 Karsten Ahnert
11 Copyright 2012 Christoph Koke
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)
19 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_DENSE_OUT_HPP_INCLUDED
20 #define BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_DENSE_OUT_HPP_INCLUDED
27 #include <boost/config.hpp> // for min/max guidelines
29 #include <boost/numeric/odeint/util/bind.hpp>
31 #include <boost/math/special_functions/binomial.hpp>
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>
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>
46 #include <boost/type_traits.hpp>
55 class Value = double ,
58 class Algebra = typename algebra_dispatcher< State >::algebra_type ,
59 class Operations = typename operations_dispatcher< State >::operations_type ,
60 class Resizer = initially_resizer
62 class bulirsch_stoer_dense_out {
67 typedef State state_type;
68 typedef Value value_type;
69 typedef Deriv deriv_type;
70 typedef Time time_type;
71 typedef Algebra algebra_type;
72 typedef Operations operations_type;
73 typedef Resizer resizer_type;
74 typedef dense_output_stepper_tag stepper_category;
76 typedef state_wrapper< state_type > wrapped_state_type;
77 typedef state_wrapper< deriv_type > wrapped_deriv_type;
79 typedef bulirsch_stoer_dense_out< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type;
81 typedef typename inverse_time< time_type >::type inv_time_type;
83 typedef std::vector< value_type > value_vector;
84 typedef std::vector< time_type > time_vector;
85 typedef std::vector< inv_time_type > inv_time_vector; //should be 1/time_type for boost.units
86 typedef std::vector< value_vector > value_matrix;
87 typedef std::vector< size_t > int_vector;
88 typedef std::vector< wrapped_state_type > state_vector_type;
89 typedef std::vector< wrapped_deriv_type > deriv_vector_type;
90 typedef std::vector< deriv_vector_type > deriv_table_type;
93 const static size_t m_k_max = 8;
97 bulirsch_stoer_dense_out(
98 value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
99 value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 ,
100 bool control_interpolation = false )
101 : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) ,
102 m_control_interpolation( control_interpolation) ,
103 m_last_step_rejected( false ) , m_first( true ) ,
104 m_current_state_x1( true ) ,
106 m_interval_sequence( m_k_max+1 ) ,
107 m_coeff( m_k_max+1 ) ,
108 m_cost( m_k_max+1 ) ,
110 m_mp_states( m_k_max+1 ) ,
111 m_derivs( m_k_max+1 ) ,
112 m_diffs( 2*m_k_max+1 ) ,
113 STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 )
115 BOOST_USING_STD_MIN();
116 BOOST_USING_STD_MAX();
118 for( unsigned short i = 0; i < m_k_max+1; i++ )
120 /* only this specific sequence allows for dense output */
121 m_interval_sequence[i] = 2 + 4*i; // 2 6 10 14 ...
122 m_derivs[i].resize( m_interval_sequence[i] );
124 m_cost[i] = m_interval_sequence[i];
126 m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
127 m_coeff[i].resize(i);
128 for( size_t k = 0 ; k < i ; ++k )
130 const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
131 m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
133 // crude estimate of optimal order
136 /* no calculation because log10 might not exist for value_type!
137 const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >( 1.0E-12 ) ) ) * 0.6 + 0.5 );
138 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 ) ));
142 for( int i = 2*(m_k_max) ; i >=0 ; i-- )
144 m_diffs[i].resize( num );
149 template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
150 controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , DerivOut &dxdt_new , time_type &dt )
152 BOOST_USING_STD_MIN();
153 BOOST_USING_STD_MAX();
156 static const value_type val1( 1.0 );
158 typename odeint::unwrap_reference< System >::type &sys = system;
162 time_vector h_opt( m_k_max+1 );
163 inv_time_vector work( m_k_max+1 );
166 time_type new_h = dt;
168 //std::cout << "t=" << t <<", dt=" << dt << ", k_opt=" << m_current_k_opt << ", first: " << m_first << std::endl;
170 for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ )
172 m_midpoint.set_steps( m_interval_sequence[k] );
175 m_midpoint.do_step( system , in , dxdt , t , out , dt , m_mp_states[k].m_v , m_derivs[k]);
179 m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt , m_mp_states[k].m_v , m_derivs[k] );
180 extrapolate( k , m_table , m_coeff , out );
181 // get error estimate
182 m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v ,
183 typename operations_type::template scale_sum2< value_type , value_type >( val1 , -val1 ) );
184 const value_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt );
185 h_opt[k] = calc_h_opt( dt , error , k );
186 work[k] = static_cast<value_type>( m_cost[k] ) / h_opt[k];
190 if( (k == m_current_k_opt-1) || m_first )
191 { // convergence before k_opt ?
196 if( (work[k] < KFAC2*work[k-1]) || (m_current_k_opt <= 2) )
198 // leave order as is (except we were in first round)
199 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 ) );
200 new_h = h_opt[k] * static_cast<value_type>( m_cost[k+1] ) / static_cast<value_type>( m_cost[k] );
202 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) ) );
207 else if( should_reject( error , k ) && !m_first )
214 if( k == m_current_k_opt )
215 { // convergence at k_opt ?
220 if( (work[k-1] < KFAC2*work[k]) )
222 m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
223 new_h = h_opt[m_current_k_opt];
225 else if( (work[k] < KFAC2*work[k-1]) && !m_last_step_rejected )
227 m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , static_cast<int>(m_current_k_opt)+1 );
228 new_h = h_opt[k]*static_cast<value_type>( m_cost[m_current_k_opt] ) / static_cast<value_type>( m_cost[k] );
230 new_h = h_opt[m_current_k_opt];
233 else if( should_reject( error , k ) )
236 new_h = h_opt[m_current_k_opt];
240 if( k == m_current_k_opt+1 )
241 { // convergence at k_opt+1 ?
245 if( work[k-2] < KFAC2*work[k-1] )
246 m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
247 if( (work[k] < KFAC2*work[m_current_k_opt]) && !m_last_step_rejected )
248 m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , static_cast<int>(k) );
249 new_h = h_opt[m_current_k_opt];
253 new_h = h_opt[m_current_k_opt];
263 //calculate dxdt for next step and dense output
264 sys( out , dxdt_new , t+dt );
266 //prepare dense output
267 value_type error = prepare_dense_output( m_k_final , in , dxdt , out , dxdt_new , dt );
269 if( error > static_cast<value_type>(10) ) // we are not as accurate for interpolation as for the steps
272 new_h = dt * pow BOOST_PREVENT_MACRO_SUBSTITUTION( error , static_cast<value_type>(-1)/(2*m_k_final+2) );
278 if( !m_last_step_rejected || (new_h < dt) )
281 m_last_step_rejected = reject;
288 template< class StateType >
289 void initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 )
291 m_resizer.adjust_size( x0 , detail::bind( &controlled_error_bs_type::template resize_impl< StateType > , detail::ref( *this ) , detail::_1 ) );
292 boost::numeric::odeint::copy( x0 , get_current_state() );
299 /* =======================================================
300 * the actual step method that should be called from outside (maybe make try_step private?)
302 template< class System >
303 std::pair< time_type , time_type > do_step( System system )
305 const size_t max_count = 1000;
309 typename odeint::unwrap_reference< System >::type &sys = system;
310 sys( get_current_state() , get_current_deriv() , m_t );
313 controlled_step_result res = fail;
318 res = try_step( system , get_current_state() , get_current_deriv() , m_t , get_old_state() , get_old_deriv() , m_dt );
320 if( count++ == max_count )
321 throw std::overflow_error( "bulirsch_stoer : too much iterations!");
323 toggle_current_state();
324 return std::make_pair( m_t_last , m_t );
327 /* performs the interpolation from a calculated step */
328 template< class StateOut >
329 void calc_state( time_type t , StateOut &x ) const
331 do_interpolation( t , x );
334 const state_type& current_state( void ) const
336 return get_current_state();
339 time_type current_time( void ) const
344 const state_type& previous_state( void ) const
346 return get_old_state();
349 time_type previous_time( void ) const
354 time_type current_time_step( void ) const
359 /** \brief Resets the internal state of the stepper. */
363 m_last_step_rejected = false;
366 template< class StateIn >
367 void adjust_size( const StateIn &x )
370 m_midpoint.adjust_size();
376 template< class StateInOut , class StateVector >
377 void extrapolate( size_t k , StateVector &table , const value_matrix &coeff , StateInOut &xest , size_t order_start_index = 0 )
378 //polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
380 static const value_type val1( 1.0 );
381 for( int j=k-1 ; j>0 ; --j )
383 m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
384 typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][j + order_start_index] ,
385 -coeff[k + order_start_index][j + order_start_index] ) );
387 m_algebra.for_each3( xest , table[0].m_v , xest ,
388 typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][0 + order_start_index] ,
389 -coeff[k + order_start_index][0 + order_start_index]) );
393 template< class StateVector >
394 void extrapolate_dense_out( size_t k , StateVector &table , const value_matrix &coeff , size_t order_start_index = 0 )
395 //polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
397 // result is written into table[0]
398 static const value_type val1( 1.0 );
399 for( int j=k ; j>1 ; --j )
401 m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
402 typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][j + order_start_index - 1] ,
403 -coeff[k + order_start_index][j + order_start_index - 1] ) );
405 m_algebra.for_each3( table[0].m_v , table[1].m_v , table[0].m_v ,
406 typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][order_start_index] ,
407 -coeff[k + order_start_index][order_start_index]) );
410 time_type calc_h_opt( time_type h , value_type error , size_t k ) const
412 BOOST_USING_STD_MIN();
413 BOOST_USING_STD_MAX();
416 value_type expo=1.0/(m_interval_sequence[k-1]);
417 value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
423 fac = STEPFAC2 / pow BOOST_PREVENT_MACRO_SUBSTITUTION( error / STEPFAC1 , expo );
424 fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( facmin/STEPFAC4 , min BOOST_PREVENT_MACRO_SUBSTITUTION( 1.0/facmin , fac ) );
429 bool in_convergence_window( size_t k ) const
431 if( (k == m_current_k_opt-1) && !m_last_step_rejected )
432 return true; // decrease order only if last step was not rejected
433 return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) );
436 bool should_reject( value_type error , size_t k ) const
438 if( k == m_current_k_opt-1 )
440 const value_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] /
441 (m_interval_sequence[0]*m_interval_sequence[0]);
442 //step will fail, criterion 17.3.17 in NR
443 return ( error > d*d );
445 else if( k == m_current_k_opt )
447 const value_type d = m_interval_sequence[m_current_k_opt+1] / m_interval_sequence[0];
448 return ( error > d*d );
453 template< class StateIn1 , class DerivIn1 , class StateIn2 , class DerivIn2 >
454 value_type prepare_dense_output( int k , const StateIn1 &x_start , const DerivIn1 &dxdt_start ,
455 const StateIn2 & /* x_end */ , const DerivIn2 & /*dxdt_end */ , time_type dt )
456 /* k is the order to which the result was approximated */
459 /* compute the coefficients of the interpolation polynomial
460 * we parametrize the interval t .. t+dt by theta = -1 .. 1
461 * we use 2k+3 values at the interval center theta=0 to obtain the interpolation coefficients
462 * the values are x(t+dt/2) and the derivatives dx/dt , ... d^(2k+2) x / dt^(2k+2) at the midpoints
463 * the derivatives are approximated via finite differences
464 * all values are obtained from interpolation of the results from the increasing orders of the midpoint calls
467 // calculate finite difference approximations to derivatives at the midpoint
468 for( int j = 0 ; j<=k ; j++ )
470 /* not working with boost units... */
471 const value_type d = m_interval_sequence[j] / ( static_cast<value_type>(2) * dt );
472 value_type f = 1.0; //factor 1/2 here because our interpolation interval has length 2 !!!
473 for( int kappa = 0 ; kappa <= 2*j+1 ; ++kappa )
475 calculate_finite_difference( j , kappa , f , dxdt_start );
480 extrapolate_dense_out( j , m_mp_states , m_coeff );
485 // extrapolate finite differences
486 for( int kappa = 0 ; kappa<=2*k+1 ; kappa++ )
488 for( int j=1 ; j<=(k-kappa/2) ; ++j )
489 extrapolate_dense_out( j , m_diffs[kappa] , m_coeff , kappa/2 );
491 // extrapolation results are now stored in m_diffs[kappa][0]
493 // divide kappa-th derivative by kappa because we need these terms for dense output interpolation
494 m_algebra.for_each1( m_diffs[kappa][0].m_v , typename operations_type::template scale< time_type >( static_cast<time_type>(d) ) );
496 d *= dt/(2*(kappa+2));
499 // 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]
501 // the error is just the highest order coefficient of the interpolation polynomial
502 // this is because we use only the midpoint theta=0 as support for the interpolation (remember that theta = -1 .. 1)
504 value_type error = 0.0;
505 if( m_control_interpolation )
507 boost::numeric::odeint::copy( m_diffs[2*k+1][0].m_v , m_err.m_v );
508 error = m_error_checker.error( m_algebra , x_start , dxdt_start , m_err.m_v , dt );
514 template< class DerivIn >
515 void calculate_finite_difference( size_t j , size_t kappa , value_type fac , const DerivIn &dxdt )
517 const int m = m_interval_sequence[j]/2-1;
518 if( kappa == 0) // no calculation required for 0th derivative of f
520 m_algebra.for_each2( m_diffs[0][j].m_v , m_derivs[j][m].m_v ,
521 typename operations_type::template scale_sum1< value_type >( fac ) );
525 // calculate the index of m_diffs for this kappa-j-combination
526 const int j_diffs = j - kappa/2;
528 m_algebra.for_each2( m_diffs[kappa][j_diffs].m_v , m_derivs[j][m+kappa].m_v ,
529 typename operations_type::template scale_sum1< value_type >( fac ) );
530 value_type sign = -1.0;
532 //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
533 for( int i = m+static_cast<int>(kappa)-2 ; i >= m-static_cast<int>(kappa) ; i -= 2 )
537 m_algebra.for_each3( m_diffs[kappa][j_diffs].m_v , m_diffs[kappa][j_diffs].m_v , m_derivs[j][i].m_v ,
538 typename operations_type::template scale_sum2< value_type , value_type >( 1.0 ,
539 sign * fac * boost::math::binomial_coefficient< value_type >( kappa , c ) ) );
543 m_algebra.for_each3( m_diffs[kappa][j_diffs].m_v , m_diffs[kappa][j_diffs].m_v , dxdt ,
544 typename operations_type::template scale_sum2< value_type , value_type >( 1.0 , sign * fac ) );
552 template< class StateOut >
553 void do_interpolation( time_type t , StateOut &out ) const
555 // interpolation polynomial is defined for theta = -1 ... 1
556 // m_k_final is the number of order-iterations done for the last step - it governs the order of the interpolation polynomial
557 const value_type theta = 2 * get_unit_value( (t - m_t_last) / (m_t - m_t_last) ) - 1;
558 // we use only values at interval center, that is theta=0, for interpolation
559 // our interpolation polynomial is thus of order 2k+2, hence we have 2k+3 terms
561 boost::numeric::odeint::copy( m_mp_states[0].m_v , out );
562 // add remaining terms: x += a_1 theta + a2 theta^2 + ... + a_{2k} theta^{2k}
563 value_type theta_pow( theta );
564 for( size_t i=0 ; i<=2*m_k_final+1 ; ++i )
566 m_algebra.for_each3( out , out , m_diffs[i][0].m_v ,
567 typename operations_type::template scale_sum2< value_type >( static_cast<value_type>(1) , theta_pow ) );
572 /* Resizer methods */
573 template< class StateIn >
574 bool resize_impl( const StateIn &x )
576 bool resized( false );
578 resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() );
579 resized |= adjust_size_by_resizeability( m_x2 , x , typename is_resizeable<state_type>::type() );
580 resized |= adjust_size_by_resizeability( m_dxdt1 , x , typename is_resizeable<state_type>::type() );
581 resized |= adjust_size_by_resizeability( m_dxdt2 , x , typename is_resizeable<state_type>::type() );
582 resized |= adjust_size_by_resizeability( m_err , x , typename is_resizeable<state_type>::type() );
584 for( size_t i = 0 ; i < m_k_max ; ++i )
585 resized |= adjust_size_by_resizeability( m_table[i] , x , typename is_resizeable<state_type>::type() );
586 for( size_t i = 0 ; i < m_k_max+1 ; ++i )
587 resized |= adjust_size_by_resizeability( m_mp_states[i] , x , typename is_resizeable<state_type>::type() );
588 for( size_t i = 0 ; i < m_k_max+1 ; ++i )
589 for( size_t j = 0 ; j < m_derivs[i].size() ; ++j )
590 resized |= adjust_size_by_resizeability( m_derivs[i][j] , x , typename is_resizeable<deriv_type>::type() );
591 for( size_t i = 0 ; i < 2*m_k_max+1 ; ++i )
592 for( size_t j = 0 ; j < m_diffs[i].size() ; ++j )
593 resized |= adjust_size_by_resizeability( m_diffs[i][j] , x , typename is_resizeable<deriv_type>::type() );
599 state_type& get_current_state( void )
601 return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
604 const state_type& get_current_state( void ) const
606 return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
609 state_type& get_old_state( void )
611 return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
614 const state_type& get_old_state( void ) const
616 return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
619 deriv_type& get_current_deriv( void )
621 return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
624 const deriv_type& get_current_deriv( void ) const
626 return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
629 deriv_type& get_old_deriv( void )
631 return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
634 const deriv_type& get_old_deriv( void ) const
636 return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
640 void toggle_current_state( void )
642 m_current_state_x1 = ! m_current_state_x1;
647 default_error_checker< value_type, algebra_type , operations_type > m_error_checker;
648 modified_midpoint_dense_out< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint;
650 bool m_control_interpolation;
652 bool m_last_step_rejected;
660 size_t m_current_k_opt;
663 algebra_type m_algebra;
665 resizer_type m_resizer;
667 wrapped_state_type m_x1 , m_x2;
668 wrapped_deriv_type m_dxdt1 , m_dxdt2;
669 wrapped_state_type m_err;
670 bool m_current_state_x1;
674 value_vector m_error; // errors of repeated midpoint steps and extrapolations
675 int_vector m_interval_sequence; // stores the successive interval counts
676 value_matrix m_coeff;
677 int_vector m_cost; // costs for interval count
679 state_vector_type m_table; // sequence of states for extrapolation
682 state_vector_type m_mp_states; // sequence of approximations of x at distance center
683 deriv_table_type m_derivs; // table of function values
684 deriv_table_type m_diffs; // table of function values
686 //wrapped_state_type m_a1 , m_a2 , m_a3 , m_a4;
688 const value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
693 /********** DOXYGEN **********/
696 * \class bulirsch_stoer_dense_out
697 * \brief The Bulirsch-Stoer algorithm.
699 * The Bulirsch-Stoer is a controlled stepper that adjusts both step size
700 * and order of the method. The algorithm uses the modified midpoint and
701 * a polynomial extrapolation compute the solution. This class also provides
702 * dense output facility.
704 * \tparam State The state type.
705 * \tparam Value The value type.
706 * \tparam Deriv The type representing the time derivative of the state.
707 * \tparam Time The time representing the independent variable - the time.
708 * \tparam Algebra The algebra type.
709 * \tparam Operations The operations type.
710 * \tparam Resizer The resizer policy type.
714 * \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 )
715 * \brief Constructs the bulirsch_stoer class, including initialization of
718 * \param eps_abs Absolute tolerance level.
719 * \param eps_rel Relative tolerance level.
720 * \param factor_x Factor for the weight of the state.
721 * \param factor_dxdt Factor for the weight of the derivative.
722 * \param control_interpolation Set true to additionally control the error of
727 * \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 )
728 * \brief Tries to perform one step.
730 * This method tries to do one step with step size dt. If the error estimate
731 * is to large, the step is rejected and the method returns fail and the
732 * step size dt is reduced. If the error estimate is acceptably small, the
733 * step is performed, success is returned and dt might be increased to make
734 * the steps as large as possible. This method also updates t if a step is
735 * performed. Also, the internal order of the stepper is adjusted if required.
737 * \param system The system function to solve, hence the r.h.s. of the ODE.
738 * It must fulfill the Simple System concept.
739 * \param in The state of the ODE which should be solved.
740 * \param dxdt The derivative of state.
741 * \param t The value of the time. Updated if the step is successful.
742 * \param out Used to store the result of the step.
743 * \param dt The step size. Updated.
744 * \return success if the step was accepted, fail otherwise.
748 * \fn bulirsch_stoer_dense_out::initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 )
749 * \brief Initializes the dense output stepper.
751 * \param x0 The initial state.
752 * \param t0 The initial time.
753 * \param dt0 The initial time step.
757 * \fn bulirsch_stoer_dense_out::do_step( System system )
758 * \brief Does one time step. This is the main method that should be used to
759 * integrate an ODE with this stepper.
760 * \note initialize has to be called before using this method to set the
761 * initial conditions x,t and the stepsize.
762 * \param system The system function to solve, hence the r.h.s. of the
763 * ordinary differential equation. It must fulfill the Simple System concept.
764 * \return Pair with start and end time of the integration step.
768 * \fn bulirsch_stoer_dense_out::calc_state( time_type t , StateOut &x ) const
769 * \brief Calculates the solution at an intermediate point within the last step
770 * \param t The time at which the solution should be calculated, has to be
771 * in the current time interval.
772 * \param x The output variable where the result is written into.
776 * \fn bulirsch_stoer_dense_out::current_state( void ) const
777 * \brief Returns the current state of the solution.
778 * \return The current state of the solution x(t).
782 * \fn bulirsch_stoer_dense_out::current_time( void ) const
783 * \brief Returns the current time of the solution.
784 * \return The current time of the solution t.
788 * \fn bulirsch_stoer_dense_out::previous_state( void ) const
789 * \brief Returns the last state of the solution.
790 * \return The last state of the solution x(t-dt).
794 * \fn bulirsch_stoer_dense_out::previous_time( void ) const
795 * \brief Returns the last time of the solution.
796 * \return The last time of the solution t-dt.
800 * \fn bulirsch_stoer_dense_out::current_time_step( void ) const
801 * \brief Returns the current step size.
802 * \return The current step size.
806 * \fn bulirsch_stoer_dense_out::adjust_size( const StateIn &x )
807 * \brief Adjust the size of all temporaries in the stepper manually.
808 * \param x A state from which the size of the temporaries to be resized is deduced.
815 #endif // BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED