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