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