1 /* [auto_generated] 2 boost/numeric/odeint/stepper/controlled_runge_kutta.hpp 3 4 [begin_description] 5 The default controlled stepper which can be used with all explicit Runge-Kutta error steppers. 6 [end_description] 7 8 Copyright 2010-2013 Karsten Ahnert 9 Copyright 2010-2013 Mario Mulansky 10 Copyright 2012 Christoph Koke 11 12 Distributed under the Boost Software License, Version 1.0. 13 (See accompanying file LICENSE_1_0.txt or 14 copy at http://www.boost.org/LICENSE_1_0.txt) 15 */ 16 17 18 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED 19 #define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED 20 21 22 23 #include <cmath> 24 25 #include <boost/config.hpp> 26 #include <boost/utility/enable_if.hpp> 27 #include <boost/type_traits/is_same.hpp> 28 29 #include <boost/numeric/odeint/util/bind.hpp> 30 #include <boost/numeric/odeint/util/unwrap_reference.hpp> 31 #include <boost/numeric/odeint/util/copy.hpp> 32 33 #include <boost/numeric/odeint/util/state_wrapper.hpp> 34 #include <boost/numeric/odeint/util/is_resizeable.hpp> 35 #include <boost/numeric/odeint/util/resizer.hpp> 36 37 #include <boost/numeric/odeint/algebra/range_algebra.hpp> 38 #include <boost/numeric/odeint/algebra/default_operations.hpp> 39 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> 40 41 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp> 42 #include <boost/numeric/odeint/stepper/stepper_categories.hpp> 43 44 namespace boost { 45 namespace numeric { 46 namespace odeint { 47 48 49 template 50 < 51 class Value , 52 class Algebra , 53 class Operations 54 > 55 class default_error_checker 56 { 57 public: 58 59 typedef Value value_type; 60 typedef Algebra algebra_type; 61 typedef Operations operations_type; 62 default_error_checker(value_type eps_abs=static_cast<value_type> (1.0e-6),value_type eps_rel=static_cast<value_type> (1.0e-6),value_type a_x=static_cast<value_type> (1),value_type a_dxdt=static_cast<value_type> (1))63 default_error_checker( 64 value_type eps_abs = static_cast< value_type >( 1.0e-6 ) , 65 value_type eps_rel = static_cast< value_type >( 1.0e-6 ) , 66 value_type a_x = static_cast< value_type >( 1 ) , 67 value_type a_dxdt = static_cast< value_type >( 1 ) ) 68 : m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt ) 69 { } 70 71 72 template< class State , class Deriv , class Err , class Time > error(const State & x_old,const Deriv & dxdt_old,Err & x_err,Time dt) const73 value_type error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const 74 { 75 return error( algebra_type() , x_old , dxdt_old , x_err , dt ); 76 } 77 78 template< class State , class Deriv , class Err , class Time > error(algebra_type & algebra,const State & x_old,const Deriv & dxdt_old,Err & x_err,Time dt) const79 value_type error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const 80 { 81 // this overwrites x_err ! 82 algebra.for_each3( x_err , x_old , dxdt_old , 83 typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * get_unit_value( dt ) ) ); 84 85 // value_type res = algebra.reduce( x_err , 86 // typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0 ) ); 87 return algebra.norm_inf( x_err ); 88 } 89 90 private: 91 92 value_type m_eps_abs; 93 value_type m_eps_rel; 94 value_type m_a_x; 95 value_type m_a_dxdt; 96 97 }; 98 99 100 101 102 103 104 105 106 /* 107 * error stepper category dispatcher 108 */ 109 template< 110 class ErrorStepper , 111 class ErrorChecker = default_error_checker< typename ErrorStepper::value_type , 112 typename ErrorStepper::algebra_type , 113 typename ErrorStepper::operations_type > , 114 class Resizer = typename ErrorStepper::resizer_type , 115 class ErrorStepperCategory = typename ErrorStepper::stepper_category 116 > 117 class controlled_runge_kutta ; 118 119 120 121 /* 122 * explicit stepper version 123 * 124 * this class introduces the following try_step overloads 125 * try_step( sys , x , t , dt ) 126 * try_step( sys , x , dxdt , t , dt ) 127 * try_step( sys , in , t , out , dt ) 128 * try_step( sys , in , dxdt , t , out , dt ) 129 */ 130 /** 131 * \brief Implements step size control for Runge-Kutta steppers with error 132 * estimation. 133 * 134 * This class implements the step size control for standard Runge-Kutta 135 * steppers with error estimation. 136 * 137 * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept. 138 * \tparam ErrorChecker The error checker 139 * \tparam Resizer The resizer policy type. 140 */ 141 template< 142 class ErrorStepper , 143 class ErrorChecker , 144 class Resizer 145 > 146 class controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > 147 { 148 149 public: 150 151 typedef ErrorStepper stepper_type; 152 typedef typename stepper_type::state_type state_type; 153 typedef typename stepper_type::value_type value_type; 154 typedef typename stepper_type::deriv_type deriv_type; 155 typedef typename stepper_type::time_type time_type; 156 typedef typename stepper_type::algebra_type algebra_type; 157 typedef typename stepper_type::operations_type operations_type; 158 typedef Resizer resizer_type; 159 typedef ErrorChecker error_checker_type; 160 typedef explicit_controlled_stepper_tag stepper_category; 161 162 #ifndef DOXYGEN_SKIP 163 typedef typename stepper_type::wrapped_state_type wrapped_state_type; 164 typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type; 165 166 typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > controlled_stepper_type; 167 #endif //DOXYGEN_SKIP 168 169 170 /** 171 * \brief Constructs the controlled Runge-Kutta stepper. 172 * \param error_checker An instance of the error checker. 173 * \param stepper An instance of the underlying stepper. 174 */ controlled_runge_kutta(const error_checker_type & error_checker=error_checker_type (),const stepper_type & stepper=stepper_type ())175 controlled_runge_kutta( 176 const error_checker_type &error_checker = error_checker_type( ) , 177 const stepper_type &stepper = stepper_type( ) 178 ) 179 : m_stepper( stepper ) , m_error_checker( error_checker ) 180 { } 181 182 183 184 /* 185 * Version 1 : try_step( sys , x , t , dt ) 186 * 187 * The overloads are needed to solve the forwarding problem 188 */ 189 /** 190 * \brief Tries to perform one step. 191 * 192 * This method tries to do one step with step size dt. If the error estimate 193 * is to large, the step is rejected and the method returns fail and the 194 * step size dt is reduced. If the error estimate is acceptably small, the 195 * step is performed, success is returned and dt might be increased to make 196 * the steps as large as possible. This method also updates t if a step is 197 * performed. 198 * 199 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 200 * Simple System concept. 201 * \param x The state of the ODE which should be solved. Overwritten if 202 * the step is successful. 203 * \param t The value of the time. Updated if the step is successful. 204 * \param dt The step size. Updated. 205 * \return success if the step was accepted, fail otherwise. 206 */ 207 template< class System , class StateInOut > try_step(System system,StateInOut & x,time_type & t,time_type & dt)208 controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt ) 209 { 210 return try_step_v1( system , x , t, dt ); 211 } 212 213 /** 214 * \brief Tries to perform one step. Solves the forwarding problem and 215 * allows for using boost range as state_type. 216 * 217 * This method tries to do one step with step size dt. If the error estimate 218 * is to large, the step is rejected and the method returns fail and the 219 * step size dt is reduced. If the error estimate is acceptably small, the 220 * step is performed, success is returned and dt might be increased to make 221 * the steps as large as possible. This method also updates t if a step is 222 * performed. 223 * 224 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 225 * Simple System concept. 226 * \param x The state of the ODE which should be solved. Overwritten if 227 * the step is successful. Can be a boost range. 228 * \param t The value of the time. Updated if the step is successful. 229 * \param dt The step size. Updated. 230 * \return success if the step was accepted, fail otherwise. 231 */ 232 template< class System , class StateInOut > try_step(System system,const StateInOut & x,time_type & t,time_type & dt)233 controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt ) 234 { 235 return try_step_v1( system , x , t, dt ); 236 } 237 238 239 240 /* 241 * Version 2 : try_step( sys , x , dxdt , t , dt ) 242 * 243 * this version does not solve the forwarding problem, boost.range can not be used 244 */ 245 /** 246 * \brief Tries to perform one step. 247 * 248 * This method tries to do one step with step size dt. If the error estimate 249 * is to large, the step is rejected and the method returns fail and the 250 * step size dt is reduced. If the error estimate is acceptably small, the 251 * step is performed, success is returned and dt might be increased to make 252 * the steps as large as possible. This method also updates t if a step is 253 * performed. 254 * 255 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 256 * Simple System concept. 257 * \param x The state of the ODE which should be solved. Overwritten if 258 * the step is successful. 259 * \param dxdt The derivative of state. 260 * \param t The value of the time. Updated if the step is successful. 261 * \param dt The step size. Updated. 262 * \return success if the step was accepted, fail otherwise. 263 */ 264 template< class System , class StateInOut , class DerivIn > try_step(System system,StateInOut & x,const DerivIn & dxdt,time_type & t,time_type & dt)265 controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt ) 266 { 267 m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); 268 controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt ); 269 if( res == success ) 270 { 271 boost::numeric::odeint::copy( m_xnew.m_v , x ); 272 } 273 return res; 274 } 275 276 /* 277 * Version 3 : try_step( sys , in , t , out , dt ) 278 * 279 * this version does not solve the forwarding problem, boost.range can not be used 280 * 281 * the disable is needed to avoid ambiguous overloads if state_type = time_type 282 */ 283 /** 284 * \brief Tries to perform one step. 285 * 286 * \note This method is disabled if state_type=time_type to avoid ambiguity. 287 * 288 * This method tries to do one step with step size dt. If the error estimate 289 * is to large, the step is rejected and the method returns fail and the 290 * step size dt is reduced. If the error estimate is acceptably small, the 291 * step is performed, success is returned and dt might be increased to make 292 * the steps as large as possible. This method also updates t if a step is 293 * performed. 294 * 295 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 296 * Simple System concept. 297 * \param in The state of the ODE which should be solved. 298 * \param t The value of the time. Updated if the step is successful. 299 * \param out Used to store the result of the step. 300 * \param dt The step size. Updated. 301 * \return success if the step was accepted, fail otherwise. 302 */ 303 template< class System , class StateIn , class StateOut > 304 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)305 try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt ) 306 { 307 typename odeint::unwrap_reference< System >::type &sys = system; 308 m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ); 309 sys( in , m_dxdt.m_v , t ); 310 return try_step( system , in , m_dxdt.m_v , t , out , dt ); 311 } 312 313 314 /* 315 * Version 4 : try_step( sys , in , dxdt , t , out , dt ) 316 * 317 * this version does not solve the forwarding problem, boost.range can not be used 318 */ 319 /** 320 * \brief Tries to perform one step. 321 * 322 * This method tries to do one step with step size dt. If the error estimate 323 * is to large, the step is rejected and the method returns fail and the 324 * step size dt is reduced. If the error estimate is acceptably small, the 325 * step is performed, success is returned and dt might be increased to make 326 * the steps as large as possible. This method also updates t if a step is 327 * performed. 328 * 329 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 330 * Simple System concept. 331 * \param in The state of the ODE which should be solved. 332 * \param dxdt The derivative of state. 333 * \param t The value of the time. Updated if the step is successful. 334 * \param out Used to store the result of the step. 335 * \param dt The step size. Updated. 336 * \return success if the step was accepted, fail otherwise. 337 */ 338 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)339 controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt ) 340 { 341 BOOST_USING_STD_MIN(); 342 BOOST_USING_STD_MAX(); 343 using std::pow; 344 345 m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ); 346 347 // do one step with error calculation 348 m_stepper.do_step( system , in , dxdt , t , out , dt , m_xerr.m_v ); 349 350 m_max_rel_error = m_error_checker.error( m_stepper.algebra() , in , dxdt , m_xerr.m_v , dt ); 351 352 if( m_max_rel_error > 1.0 ) 353 { 354 // error too large - decrease dt ,limit scaling factor to 0.2 and reset state 355 dt *= max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) * 356 pow( m_max_rel_error , static_cast<value_type>(-1) / ( m_stepper.error_order() - 1 ) ) ) , 357 static_cast<value_type>( static_cast<value_type>(1)/static_cast<value_type> (5) ) ); 358 return fail; 359 } 360 else 361 { 362 if( m_max_rel_error < 0.5 ) 363 { 364 // error should be > 0 365 m_max_rel_error = max BOOST_PREVENT_MACRO_SUBSTITUTION ( 366 static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(m_stepper.stepper_order()) ) ) , 367 m_max_rel_error ); 368 //error too small - increase dt and keep the evolution and limit scaling factor to 5.0 369 t += dt; 370 dt *= static_cast<value_type>(9)/static_cast<value_type>(10) * pow( m_max_rel_error , 371 static_cast<value_type>(-1) / m_stepper.stepper_order() ); 372 return success; 373 } 374 else 375 { 376 t += dt; 377 return success; 378 } 379 } 380 } 381 382 /** 383 * \brief Returns the error of the last step. 384 * 385 * returns The last error of the step. 386 */ last_error(void) const387 value_type last_error( void ) const 388 { 389 return m_max_rel_error; 390 } 391 392 393 /** 394 * \brief Adjust the size of all temporaries in the stepper manually. 395 * \param x A state from which the size of the temporaries to be resized is deduced. 396 */ 397 template< class StateType > adjust_size(const StateType & x)398 void adjust_size( const StateType &x ) 399 { 400 resize_m_xerr_impl( x ); 401 resize_m_dxdt_impl( x ); 402 resize_m_xnew_impl( x ); 403 m_stepper.adjust_size( x ); 404 } 405 406 /** 407 * \brief Returns the instance of the underlying stepper. 408 * \returns The instance of the underlying stepper. 409 */ stepper(void)410 stepper_type& stepper( void ) 411 { 412 return m_stepper; 413 } 414 415 /** 416 * \brief Returns the instance of the underlying stepper. 417 * \returns The instance of the underlying stepper. 418 */ stepper(void) const419 const stepper_type& stepper( void ) const 420 { 421 return m_stepper; 422 } 423 424 private: 425 426 427 template< class System , class StateInOut > try_step_v1(System system,StateInOut & x,time_type & t,time_type & dt)428 controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt ) 429 { 430 typename odeint::unwrap_reference< System >::type &sys = system; 431 m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); 432 sys( x , m_dxdt.m_v ,t ); 433 return try_step( system , x , m_dxdt.m_v , t , dt ); 434 } 435 436 template< class StateIn > resize_m_xerr_impl(const StateIn & x)437 bool resize_m_xerr_impl( const StateIn &x ) 438 { 439 return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() ); 440 } 441 442 template< class StateIn > resize_m_dxdt_impl(const StateIn & x)443 bool resize_m_dxdt_impl( const StateIn &x ) 444 { 445 return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() ); 446 } 447 448 template< class StateIn > resize_m_xnew_impl(const StateIn & x)449 bool resize_m_xnew_impl( const StateIn &x ) 450 { 451 return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() ); 452 } 453 454 455 456 stepper_type m_stepper; 457 error_checker_type m_error_checker; 458 459 resizer_type m_dxdt_resizer; 460 resizer_type m_xerr_resizer; 461 resizer_type m_xnew_resizer; 462 463 wrapped_deriv_type m_dxdt; 464 wrapped_state_type m_xerr; 465 wrapped_state_type m_xnew; 466 value_type m_max_rel_error; 467 }; 468 469 470 471 472 473 474 475 476 477 478 /* 479 * explicit stepper fsal version 480 * 481 * the class introduces the following try_step overloads 482 * try_step( sys , x , t , dt ) 483 * try_step( sys , in , t , out , dt ) 484 * try_step( sys , x , dxdt , t , dt ) 485 * try_step( sys , in , dxdt_in , t , out , dxdt_out , dt ) 486 */ 487 /** 488 * \brief Implements step size control for Runge-Kutta FSAL steppers with 489 * error estimation. 490 * 491 * This class implements the step size control for FSAL Runge-Kutta 492 * steppers with error estimation. 493 * 494 * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept. 495 * \tparam ErrorChecker The error checker 496 * \tparam Resizer The resizer policy type. 497 */ 498 template< 499 class ErrorStepper , 500 class ErrorChecker , 501 class Resizer 502 > 503 class controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_fsal_tag > 504 { 505 506 public: 507 508 typedef ErrorStepper stepper_type; 509 typedef typename stepper_type::state_type state_type; 510 typedef typename stepper_type::value_type value_type; 511 typedef typename stepper_type::deriv_type deriv_type; 512 typedef typename stepper_type::time_type time_type; 513 typedef typename stepper_type::algebra_type algebra_type; 514 typedef typename stepper_type::operations_type operations_type; 515 typedef Resizer resizer_type; 516 typedef ErrorChecker error_checker_type; 517 typedef explicit_controlled_stepper_fsal_tag stepper_category; 518 519 #ifndef DOXYGEN_SKIP 520 typedef typename stepper_type::wrapped_state_type wrapped_state_type; 521 typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type; 522 523 typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > controlled_stepper_type; 524 #endif // DOXYGEN_SKIP 525 526 /** 527 * \brief Constructs the controlled Runge-Kutta stepper. 528 * \param error_checker An instance of the error checker. 529 * \param stepper An instance of the underlying stepper. 530 */ controlled_runge_kutta(const error_checker_type & error_checker=error_checker_type (),const stepper_type & stepper=stepper_type ())531 controlled_runge_kutta( 532 const error_checker_type &error_checker = error_checker_type() , 533 const stepper_type &stepper = stepper_type() 534 ) 535 : m_stepper( stepper ) , m_error_checker( error_checker ) , 536 m_first_call( true ) 537 { } 538 539 /* 540 * Version 1 : try_step( sys , x , t , dt ) 541 * 542 * The two overloads are needed in order to solve the forwarding problem 543 */ 544 /** 545 * \brief Tries to perform one step. 546 * 547 * This method tries to do one step with step size dt. If the error estimate 548 * is to large, the step is rejected and the method returns fail and the 549 * step size dt is reduced. If the error estimate is acceptably small, the 550 * step is performed, success is returned and dt might be increased to make 551 * the steps as large as possible. This method also updates t if a step is 552 * performed. 553 * 554 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 555 * Simple System concept. 556 * \param x The state of the ODE which should be solved. Overwritten if 557 * the step is successful. 558 * \param t The value of the time. Updated if the step is successful. 559 * \param dt The step size. Updated. 560 * \return success if the step was accepted, fail otherwise. 561 */ 562 template< class System , class StateInOut > try_step(System system,StateInOut & x,time_type & t,time_type & dt)563 controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt ) 564 { 565 return try_step_v1( system , x , t , dt ); 566 } 567 568 569 /** 570 * \brief Tries to perform one step. Solves the forwarding problem and 571 * allows for using boost range as state_type. 572 * 573 * This method tries to do one step with step size dt. If the error estimate 574 * is to large, the step is rejected and the method returns fail and the 575 * step size dt is reduced. If the error estimate is acceptably small, the 576 * step is performed, success is returned and dt might be increased to make 577 * the steps as large as possible. This method also updates t if a step is 578 * performed. 579 * 580 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 581 * Simple System concept. 582 * \param x The state of the ODE which should be solved. Overwritten if 583 * the step is successful. Can be a boost range. 584 * \param t The value of the time. Updated if the step is successful. 585 * \param dt The step size. Updated. 586 * \return success if the step was accepted, fail otherwise. 587 */ 588 template< class System , class StateInOut > try_step(System system,const StateInOut & x,time_type & t,time_type & dt)589 controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt ) 590 { 591 return try_step_v1( system , x , t , dt ); 592 } 593 594 595 596 /* 597 * Version 2 : try_step( sys , in , t , out , dt ); 598 * 599 * This version does not solve the forwarding problem, boost::range can not be used. 600 * 601 * The disabler is needed to solve ambiguous overloads 602 */ 603 /** 604 * \brief Tries to perform one step. 605 * 606 * \note This method is disabled if state_type=time_type to avoid ambiguity. 607 * 608 * This method tries to do one step with step size dt. If the error estimate 609 * is to large, the step is rejected and the method returns fail and the 610 * step size dt is reduced. If the error estimate is acceptably small, the 611 * step is performed, success is returned and dt might be increased to make 612 * the steps as large as possible. This method also updates t if a step is 613 * performed. 614 * 615 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 616 * Simple System concept. 617 * \param in The state of the ODE which should be solved. 618 * \param t The value of the time. Updated if the step is successful. 619 * \param out Used to store the result of the step. 620 * \param dt The step size. Updated. 621 * \return success if the step was accepted, fail otherwise. 622 */ 623 template< class System , class StateIn , class StateOut > 624 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)625 try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt ) 626 { 627 if( m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) || m_first_call ) 628 { 629 initialize( system , in , t ); 630 } 631 return try_step( system , in , m_dxdt.m_v , t , out , dt ); 632 } 633 634 635 /* 636 * Version 3 : try_step( sys , x , dxdt , t , dt ) 637 * 638 * This version does not solve the forwarding problem, boost::range can not be used. 639 */ 640 /** 641 * \brief Tries to perform one step. 642 * 643 * This method tries to do one step with step size dt. If the error estimate 644 * is to large, the step is rejected and the method returns fail and the 645 * step size dt is reduced. If the error estimate is acceptably small, the 646 * step is performed, success is returned and dt might be increased to make 647 * the steps as large as possible. This method also updates t if a step is 648 * performed. 649 * 650 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 651 * Simple System concept. 652 * \param x The state of the ODE which should be solved. Overwritten if 653 * the step is successful. 654 * \param dxdt The derivative of state. 655 * \param t The value of the time. Updated if the step is successful. 656 * \param dt The step size. Updated. 657 * \return success if the step was accepted, fail otherwise. 658 */ 659 template< class System , class StateInOut , class DerivInOut > try_step(System system,StateInOut & x,DerivInOut & dxdt,time_type & t,time_type & dt)660 controlled_step_result try_step( System system , StateInOut &x , DerivInOut &dxdt , time_type &t , time_type &dt ) 661 { 662 m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); 663 m_dxdt_new_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_new_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); 664 controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , m_dxdtnew.m_v , dt ); 665 if( res == success ) 666 { 667 boost::numeric::odeint::copy( m_xnew.m_v , x ); 668 boost::numeric::odeint::copy( m_dxdtnew.m_v , dxdt ); 669 } 670 return res; 671 } 672 673 674 /* 675 * Version 4 : try_step( sys , in , dxdt_in , t , out , dxdt_out , dt ) 676 * 677 * This version does not solve the forwarding problem, boost::range can not be used. 678 */ 679 /** 680 * \brief Tries to perform one step. 681 * 682 * This method tries to do one step with step size dt. If the error estimate 683 * is to large, the step is rejected and the method returns fail and the 684 * step size dt is reduced. If the error estimate is acceptably small, the 685 * step is performed, success is returned and dt might be increased to make 686 * the steps as large as possible. This method also updates t if a step is 687 * performed. 688 * 689 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 690 * Simple System concept. 691 * \param in The state of the ODE which should be solved. 692 * \param dxdt The derivative of state. 693 * \param t The value of the time. Updated if the step is successful. 694 * \param out Used to store the result of the step. 695 * \param dt The step size. Updated. 696 * \return success if the step was accepted, fail otherwise. 697 */ 698 template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut > try_step(System system,const StateIn & in,const DerivIn & dxdt_in,time_type & t,StateOut & out,DerivOut & dxdt_out,time_type & dt)699 controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t , 700 StateOut &out , DerivOut &dxdt_out , time_type &dt ) 701 { 702 BOOST_USING_STD_MIN(); 703 BOOST_USING_STD_MAX(); 704 705 using std::pow; 706 707 m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ); 708 709 //fsal: m_stepper.get_dxdt( dxdt ); 710 //fsal: m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err ); 711 m_stepper.do_step( system , in , dxdt_in , t , out , dxdt_out , dt , m_xerr.m_v ); 712 713 // this potentially overwrites m_x_err! (standard_error_checker does, at least) 714 value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt_in , m_xerr.m_v , dt ); 715 716 if( max_rel_err > 1.0 ) 717 { 718 // error too large - decrease dt ,limit scaling factor to 0.2 and reset state 719 dt *= max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) * 720 pow( max_rel_err , static_cast<value_type>(-1) / ( m_stepper.error_order() - 1 ) ) ) , 721 static_cast<value_type>( static_cast<value_type>(1)/static_cast<value_type> (5)) ); 722 return fail; 723 } 724 else 725 { 726 if( max_rel_err < 0.5 ) 727 { //error too small - increase dt and keep the evolution and limit scaling factor to 5.0 728 // error should be > 0 729 max_rel_err = max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(m_stepper.stepper_order()) ) ) , 730 max_rel_err ); 731 t += dt; 732 dt *= static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) * pow( max_rel_err , static_cast<value_type>(-1) / m_stepper.stepper_order() ) ); 733 return success; 734 } 735 else 736 { 737 t += dt; 738 return success; 739 } 740 } 741 } 742 743 744 /** 745 * \brief Resets the internal state of the underlying FSAL stepper. 746 */ reset(void)747 void reset( void ) 748 { 749 m_first_call = true; 750 } 751 752 /** 753 * \brief Initializes the internal state storing an internal copy of the derivative. 754 * 755 * \param deriv The initial derivative of the ODE. 756 */ 757 template< class DerivIn > initialize(const DerivIn & deriv)758 void initialize( const DerivIn &deriv ) 759 { 760 boost::numeric::odeint::copy( deriv , m_dxdt.m_v ); 761 m_first_call = false; 762 } 763 764 /** 765 * \brief Initializes the internal state storing an internal copy of the derivative. 766 * 767 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 768 * Simple System concept. 769 * \param x The initial state of the ODE which should be solved. 770 * \param t The initial time. 771 */ 772 template< class System , class StateIn > initialize(System system,const StateIn & x,time_type t)773 void initialize( System system , const StateIn &x , time_type t ) 774 { 775 typename odeint::unwrap_reference< System >::type &sys = system; 776 sys( x , m_dxdt.m_v , t ); 777 m_first_call = false; 778 } 779 780 /** 781 * \brief Returns true if the stepper has been initialized, false otherwise. 782 * 783 * \return true, if the stepper has been initialized, false otherwise. 784 */ is_initialized(void) const785 bool is_initialized( void ) const 786 { 787 return ! m_first_call; 788 } 789 790 791 /** 792 * \brief Adjust the size of all temporaries in the stepper manually. 793 * \param x A state from which the size of the temporaries to be resized is deduced. 794 */ 795 template< class StateType > adjust_size(const StateType & x)796 void adjust_size( const StateType &x ) 797 { 798 resize_m_xerr_impl( x ); 799 resize_m_dxdt_impl( x ); 800 resize_m_dxdt_new_impl( x ); 801 resize_m_xnew_impl( x ); 802 } 803 804 805 /** 806 * \brief Returns the instance of the underlying stepper. 807 * \returns The instance of the underlying stepper. 808 */ stepper(void)809 stepper_type& stepper( void ) 810 { 811 return m_stepper; 812 } 813 814 /** 815 * \brief Returns the instance of the underlying stepper. 816 * \returns The instance of the underlying stepper. 817 */ stepper(void) const818 const stepper_type& stepper( void ) const 819 { 820 return m_stepper; 821 } 822 823 824 825 private: 826 827 828 template< class StateIn > resize_m_xerr_impl(const StateIn & x)829 bool resize_m_xerr_impl( const StateIn &x ) 830 { 831 return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() ); 832 } 833 834 template< class StateIn > resize_m_dxdt_impl(const StateIn & x)835 bool resize_m_dxdt_impl( const StateIn &x ) 836 { 837 return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() ); 838 } 839 840 template< class StateIn > resize_m_dxdt_new_impl(const StateIn & x)841 bool resize_m_dxdt_new_impl( const StateIn &x ) 842 { 843 return adjust_size_by_resizeability( m_dxdtnew , x , typename is_resizeable<deriv_type>::type() ); 844 } 845 846 template< class StateIn > resize_m_xnew_impl(const StateIn & x)847 bool resize_m_xnew_impl( const StateIn &x ) 848 { 849 return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() ); 850 } 851 852 853 template< class System , class StateInOut > try_step_v1(System system,StateInOut & x,time_type & t,time_type & dt)854 controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt ) 855 { 856 if( m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ) || m_first_call ) 857 { 858 initialize( system , x , t ); 859 } 860 return try_step( system , x , m_dxdt.m_v , t , dt ); 861 } 862 863 864 stepper_type m_stepper; 865 error_checker_type m_error_checker; 866 867 resizer_type m_dxdt_resizer; 868 resizer_type m_xerr_resizer; 869 resizer_type m_xnew_resizer; 870 resizer_type m_dxdt_new_resizer; 871 872 wrapped_deriv_type m_dxdt; 873 wrapped_state_type m_xerr; 874 wrapped_state_type m_xnew; 875 wrapped_deriv_type m_dxdtnew; 876 bool m_first_call; 877 }; 878 879 880 /********** DOXYGEN **********/ 881 882 /**** DEFAULT ERROR CHECKER ****/ 883 884 /** 885 * \class default_error_checker 886 * \brief The default error checker to be used with Runge-Kutta error steppers 887 * 888 * This class provides the default mechanism to compare the error estimates 889 * reported by Runge-Kutta error steppers with user defined error bounds. 890 * It is used by the controlled_runge_kutta steppers. 891 * 892 * \tparam Value The value type. 893 * \tparam Algebra The algebra type. 894 * \tparam Operations The operations type. 895 */ 896 897 /** 898 * \fn default_error_checker( value_type eps_abs , value_type eps_rel , value_type a_x , value_type a_dxdt ) 899 * \brief Constructs the error checker. 900 * 901 * The error is calculated as follows: ???? 902 * 903 * \param eps_abs Absolute tolerance level. 904 * \param eps_rel Relative tolerance level. 905 * \param a_x Factor for the weight of the state. 906 * \param a_dxdt Factor for the weight of the derivative. 907 */ 908 909 /** 910 * \fn error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const 911 * \brief Calculates the error level. 912 * 913 * If the returned error level is greater than 1, the estimated error was 914 * larger than the permitted error bounds and the step should be repeated 915 * with a smaller step size. 916 * 917 * \param x_old State at the beginning of the step. 918 * \param dxdt_old Derivative at the beginning of the step. 919 * \param x_err Error estimate. 920 * \param dt Time step. 921 * \return error 922 */ 923 924 /** 925 * \fn error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const 926 * \brief Calculates the error level using a given algebra. 927 * 928 * If the returned error level is greater than 1, the estimated error was 929 * larger than the permitted error bounds and the step should be repeated 930 * with a smaller step size. 931 * 932 * \param algebra The algebra used for calculation of the error. 933 * \param x_old State at the beginning of the step. 934 * \param dxdt_old Derivative at the beginning of the step. 935 * \param x_err Error estimate. 936 * \param dt Time step. 937 * \return error 938 */ 939 940 941 } // odeint 942 } // numeric 943 } // boost 944 945 946 #endif // BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED 947