2 boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
5 The default controlled stepper which can be used with all explicit Runge-Kutta error steppers.
8 Copyright 2010-2013 Karsten Ahnert
9 Copyright 2010-2015 Mario Mulansky
10 Copyright 2012 Christoph Koke
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)
18 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED
19 #define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED
25 #include <boost/config.hpp>
26 #include <boost/utility/enable_if.hpp>
27 #include <boost/type_traits/is_same.hpp>
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>
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 #include <boost/numeric/odeint/util/detail/less_with_sign.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>
42 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
43 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
56 class default_error_checker
60 typedef Value value_type;
61 typedef Algebra algebra_type;
62 typedef Operations operations_type;
64 default_error_checker(
65 value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
66 value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
67 value_type a_x = static_cast< value_type >( 1 ) ,
68 value_type a_dxdt = static_cast< value_type >( 1 ))
69 : m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
73 template< class State , class Deriv , class Err, class Time >
74 value_type error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
76 return error( algebra_type() , x_old , dxdt_old , x_err , dt );
79 template< class State , class Deriv , class Err, class Time >
80 value_type error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
83 // this overwrites x_err !
84 algebra.for_each3( x_err , x_old , dxdt_old ,
85 typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * abs(get_unit_value( dt )) ) );
87 // value_type res = algebra.reduce( x_err ,
88 // typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0 ) );
89 return algebra.norm_inf( x_err );
102 template< typename Value, typename Time >
103 class default_step_adjuster
106 typedef Time time_type;
107 typedef Value value_type;
109 default_step_adjuster(const time_type max_dt=static_cast<time_type>(0))
114 time_type decrease_step(time_type dt, const value_type error, const int error_order) const
116 // returns the decreased time step
117 BOOST_USING_STD_MIN();
118 BOOST_USING_STD_MAX();
122 BOOST_PREVENT_MACRO_SUBSTITUTION(
123 static_cast<value_type>( static_cast<value_type>(9) / static_cast<value_type>(10) *
124 pow(error, static_cast<value_type>(-1) / (error_order - 1))),
125 static_cast<value_type>( static_cast<value_type>(1) / static_cast<value_type> (5)));
126 if(m_max_dt != static_cast<time_type >(0))
127 // limit to maximal stepsize even when decreasing
128 dt = detail::min_abs(dt, m_max_dt);
132 time_type increase_step(time_type dt, value_type error, const int stepper_order) const
134 // returns the increased time step
135 BOOST_USING_STD_MIN();
136 BOOST_USING_STD_MAX();
139 // adjust the size if dt is smaller than max_dt (providede max_dt is not zero)
142 // error should be > 0
143 error = max BOOST_PREVENT_MACRO_SUBSTITUTION (
144 static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(stepper_order) ) ) ,
146 // time_type dt_old = dt; unused variable warning
147 //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
148 dt *= static_cast<value_type>(9)/static_cast<value_type>(10) *
149 pow(error, static_cast<value_type>(-1) / stepper_order);
150 if(m_max_dt != static_cast<time_type >(0))
151 // limit to maximal stepsize
152 dt = detail::min_abs(dt, m_max_dt);
157 bool check_step_size_limit(const time_type dt)
159 if(m_max_dt != static_cast<time_type >(0))
160 return detail::less_eq_with_sign(dt, m_max_dt, dt);
164 time_type get_max_dt() { return m_max_dt; }
173 * error stepper category dispatcher
177 class ErrorChecker = default_error_checker< typename ErrorStepper::value_type ,
178 typename ErrorStepper::algebra_type ,
179 typename ErrorStepper::operations_type > ,
180 class StepAdjuster = default_step_adjuster< typename ErrorStepper::value_type ,
181 typename ErrorStepper::time_type > ,
182 class Resizer = typename ErrorStepper::resizer_type ,
183 class ErrorStepperCategory = typename ErrorStepper::stepper_category
185 class controlled_runge_kutta ;
190 * explicit stepper version
192 * this class introduces the following try_step overloads
193 * try_step( sys , x , t , dt )
194 * try_step( sys , x , dxdt , t , dt )
195 * try_step( sys , in , t , out , dt )
196 * try_step( sys , in , dxdt , t , out , dt )
199 * \brief Implements step size control for Runge-Kutta steppers with error
202 * This class implements the step size control for standard Runge-Kutta
203 * steppers with error estimation.
205 * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept.
206 * \tparam ErrorChecker The error checker
207 * \tparam Resizer The resizer policy type.
215 class controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster, Resizer ,
216 explicit_error_stepper_tag >
221 typedef ErrorStepper stepper_type;
222 typedef typename stepper_type::state_type state_type;
223 typedef typename stepper_type::value_type value_type;
224 typedef typename stepper_type::deriv_type deriv_type;
225 typedef typename stepper_type::time_type time_type;
226 typedef typename stepper_type::algebra_type algebra_type;
227 typedef typename stepper_type::operations_type operations_type;
228 typedef Resizer resizer_type;
229 typedef ErrorChecker error_checker_type;
230 typedef StepAdjuster step_adjuster_type;
231 typedef explicit_controlled_stepper_tag stepper_category;
234 typedef typename stepper_type::wrapped_state_type wrapped_state_type;
235 typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
237 typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster ,
238 Resizer , explicit_error_stepper_tag > controlled_stepper_type;
239 #endif //DOXYGEN_SKIP
243 * \brief Constructs the controlled Runge-Kutta stepper.
244 * \param error_checker An instance of the error checker.
245 * \param stepper An instance of the underlying stepper.
247 controlled_runge_kutta(
248 const error_checker_type &error_checker = error_checker_type( ) ,
249 const step_adjuster_type &step_adjuster = step_adjuster_type() ,
250 const stepper_type &stepper = stepper_type( )
252 : m_stepper(stepper), m_error_checker(error_checker) , m_step_adjuster(step_adjuster)
258 * Version 1 : try_step( sys , x , t , dt )
260 * The overloads are needed to solve the forwarding problem
263 * \brief Tries to perform one step.
265 * This method tries to do one step with step size dt. If the error estimate
266 * is to large, the step is rejected and the method returns fail and the
267 * step size dt is reduced. If the error estimate is acceptably small, the
268 * step is performed, success is returned and dt might be increased to make
269 * the steps as large as possible. This method also updates t if a step is
272 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
273 * Simple System concept.
274 * \param x The state of the ODE which should be solved. Overwritten if
275 * the step is successful.
276 * \param t The value of the time. Updated if the step is successful.
277 * \param dt The step size. Updated.
278 * \return success if the step was accepted, fail otherwise.
280 template< class System , class StateInOut >
281 controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
283 return try_step_v1( system , x , t, dt );
287 * \brief Tries to perform one step. Solves the forwarding problem and
288 * allows for using boost range as state_type.
290 * This method tries to do one step with step size dt. If the error estimate
291 * is to large, the step is rejected and the method returns fail and the
292 * step size dt is reduced. If the error estimate is acceptably small, the
293 * step is performed, success is returned and dt might be increased to make
294 * the steps as large as possible. This method also updates t if a step is
297 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
298 * Simple System concept.
299 * \param x The state of the ODE which should be solved. Overwritten if
300 * the step is successful. Can be a boost range.
301 * \param t The value of the time. Updated if the step is successful.
302 * \param dt The step size. Updated.
303 * \return success if the step was accepted, fail otherwise.
305 template< class System , class StateInOut >
306 controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
308 return try_step_v1( system , x , t, dt );
314 * Version 2 : try_step( sys , x , dxdt , t , dt )
316 * this version does not solve the forwarding problem, boost.range can not be used
319 * \brief Tries to perform one step.
321 * This method tries to do one step with step size dt. If the error estimate
322 * is to large, the step is rejected and the method returns fail and the
323 * step size dt is reduced. If the error estimate is acceptably small, the
324 * step is performed, success is returned and dt might be increased to make
325 * the steps as large as possible. This method also updates t if a step is
328 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
329 * Simple System concept.
330 * \param x The state of the ODE which should be solved. Overwritten if
331 * the step is successful.
332 * \param dxdt The derivative of state.
333 * \param t The value of the time. Updated if the step is successful.
334 * \param dt The step size. Updated.
335 * \return success if the step was accepted, fail otherwise.
337 template< class System , class StateInOut , class DerivIn >
338 controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt )
340 m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
341 controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt );
344 boost::numeric::odeint::copy( m_xnew.m_v , x );
350 * Version 3 : try_step( sys , in , t , out , dt )
352 * this version does not solve the forwarding problem, boost.range can not be used
354 * the disable is needed to avoid ambiguous overloads if state_type = time_type
357 * \brief Tries to perform one step.
359 * \note This method is disabled if state_type=time_type to avoid ambiguity.
361 * This method tries to do one step with step size dt. If the error estimate
362 * is to large, the step is rejected and the method returns fail and the
363 * step size dt is reduced. If the error estimate is acceptably small, the
364 * step is performed, success is returned and dt might be increased to make
365 * the steps as large as possible. This method also updates t if a step is
368 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
369 * Simple System concept.
370 * \param in The state of the ODE which should be solved.
371 * \param t The value of the time. Updated if the step is successful.
372 * \param out Used to store the result of the step.
373 * \param dt The step size. Updated.
374 * \return success if the step was accepted, fail otherwise.
376 template< class System , class StateIn , class StateOut >
377 typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type
378 try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
380 typename odeint::unwrap_reference< System >::type &sys = system;
381 m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
382 sys( in , m_dxdt.m_v , t );
383 return try_step( system , in , m_dxdt.m_v , t , out , dt );
388 * Version 4 : try_step( sys , in , dxdt , t , out , dt )
390 * this version does not solve the forwarding problem, boost.range can not be used
393 * \brief Tries to perform one step.
395 * This method tries to do one step with step size dt. If the error estimate
396 * is to large, the step is rejected and the method returns fail and the
397 * step size dt is reduced. If the error estimate is acceptably small, the
398 * step is performed, success is returned and dt might be increased to make
399 * the steps as large as possible. This method also updates t if a step is
402 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
403 * Simple System concept.
404 * \param in The state of the ODE which should be solved.
405 * \param dxdt The derivative of state.
406 * \param t The value of the time. Updated if the step is successful.
407 * \param out Used to store the result of the step.
408 * \param dt The step size. Updated.
409 * \return success if the step was accepted, fail otherwise.
411 template< class System , class StateIn , class DerivIn , class StateOut >
412 controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
414 if( !m_step_adjuster.check_step_size_limit(dt) )
416 // given dt was above step size limit - adjust and return fail;
417 dt = m_step_adjuster.get_max_dt();
421 m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
423 // do one step with error calculation
424 m_stepper.do_step( system , in , dxdt , t , out , dt , m_xerr.m_v );
426 value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt , m_xerr.m_v , dt );
428 if( max_rel_err > 1.0 )
430 // error too big, decrease step size and reject this step
431 dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
435 // otherwise, increase step size and accept
437 dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
443 * \brief Adjust the size of all temporaries in the stepper manually.
444 * \param x A state from which the size of the temporaries to be resized is deduced.
446 template< class StateType >
447 void adjust_size( const StateType &x )
449 resize_m_xerr_impl( x );
450 resize_m_dxdt_impl( x );
451 resize_m_xnew_impl( x );
452 m_stepper.adjust_size( x );
456 * \brief Returns the instance of the underlying stepper.
457 * \returns The instance of the underlying stepper.
459 stepper_type& stepper( void )
465 * \brief Returns the instance of the underlying stepper.
466 * \returns The instance of the underlying stepper.
468 const stepper_type& stepper( void ) const
476 template< class System , class StateInOut >
477 controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
479 typename odeint::unwrap_reference< System >::type &sys = system;
480 m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
481 sys( x , m_dxdt.m_v ,t );
482 return try_step( system , x , m_dxdt.m_v , t , dt );
485 template< class StateIn >
486 bool resize_m_xerr_impl( const StateIn &x )
488 return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() );
491 template< class StateIn >
492 bool resize_m_dxdt_impl( const StateIn &x )
494 return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
497 template< class StateIn >
498 bool resize_m_xnew_impl( const StateIn &x )
500 return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
505 stepper_type m_stepper;
506 error_checker_type m_error_checker;
507 step_adjuster_type m_step_adjuster;
509 resizer_type m_dxdt_resizer;
510 resizer_type m_xerr_resizer;
511 resizer_type m_xnew_resizer;
513 wrapped_deriv_type m_dxdt;
514 wrapped_state_type m_xerr;
515 wrapped_state_type m_xnew;
528 * explicit stepper fsal version
530 * the class introduces the following try_step overloads
531 * try_step( sys , x , t , dt )
532 * try_step( sys , in , t , out , dt )
533 * try_step( sys , x , dxdt , t , dt )
534 * try_step( sys , in , dxdt_in , t , out , dxdt_out , dt )
537 * \brief Implements step size control for Runge-Kutta FSAL steppers with
540 * This class implements the step size control for FSAL Runge-Kutta
541 * steppers with error estimation.
543 * \tparam ErrorStepper The stepper type with error estimation, has to fulfill the ErrorStepper concept.
544 * \tparam ErrorChecker The error checker
545 * \tparam Resizer The resizer policy type.
553 class controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster , Resizer , explicit_error_stepper_fsal_tag >
558 typedef ErrorStepper stepper_type;
559 typedef typename stepper_type::state_type state_type;
560 typedef typename stepper_type::value_type value_type;
561 typedef typename stepper_type::deriv_type deriv_type;
562 typedef typename stepper_type::time_type time_type;
563 typedef typename stepper_type::algebra_type algebra_type;
564 typedef typename stepper_type::operations_type operations_type;
565 typedef Resizer resizer_type;
566 typedef ErrorChecker error_checker_type;
567 typedef StepAdjuster step_adjuster_type;
568 typedef explicit_controlled_stepper_fsal_tag stepper_category;
571 typedef typename stepper_type::wrapped_state_type wrapped_state_type;
572 typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
574 typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
575 #endif // DOXYGEN_SKIP
578 * \brief Constructs the controlled Runge-Kutta stepper.
579 * \param error_checker An instance of the error checker.
580 * \param stepper An instance of the underlying stepper.
582 controlled_runge_kutta(
583 const error_checker_type &error_checker = error_checker_type() ,
584 const step_adjuster_type &step_adjuster = step_adjuster_type() ,
585 const stepper_type &stepper = stepper_type()
587 : m_stepper( stepper ) , m_error_checker( error_checker ) , m_step_adjuster(step_adjuster) ,
592 * Version 1 : try_step( sys , x , t , dt )
594 * The two overloads are needed in order to solve the forwarding problem
597 * \brief Tries to perform one step.
599 * This method tries to do one step with step size dt. If the error estimate
600 * is to large, the step is rejected and the method returns fail and the
601 * step size dt is reduced. If the error estimate is acceptably small, the
602 * step is performed, success is returned and dt might be increased to make
603 * the steps as large as possible. This method also updates t if a step is
606 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
607 * Simple System concept.
608 * \param x The state of the ODE which should be solved. Overwritten if
609 * the step is successful.
610 * \param t The value of the time. Updated if the step is successful.
611 * \param dt The step size. Updated.
612 * \return success if the step was accepted, fail otherwise.
614 template< class System , class StateInOut >
615 controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
617 return try_step_v1( system , x , t , dt );
622 * \brief Tries to perform one step. Solves the forwarding problem and
623 * allows for using boost range as state_type.
625 * This method tries to do one step with step size dt. If the error estimate
626 * is to large, the step is rejected and the method returns fail and the
627 * step size dt is reduced. If the error estimate is acceptably small, the
628 * step is performed, success is returned and dt might be increased to make
629 * the steps as large as possible. This method also updates t if a step is
632 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
633 * Simple System concept.
634 * \param x The state of the ODE which should be solved. Overwritten if
635 * the step is successful. Can be a boost range.
636 * \param t The value of the time. Updated if the step is successful.
637 * \param dt The step size. Updated.
638 * \return success if the step was accepted, fail otherwise.
640 template< class System , class StateInOut >
641 controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
643 return try_step_v1( system , x , t , dt );
649 * Version 2 : try_step( sys , in , t , out , dt );
651 * This version does not solve the forwarding problem, boost::range can not be used.
653 * The disabler is needed to solve ambiguous overloads
656 * \brief Tries to perform one step.
658 * \note This method is disabled if state_type=time_type to avoid ambiguity.
660 * This method tries to do one step with step size dt. If the error estimate
661 * is to large, the step is rejected and the method returns fail and the
662 * step size dt is reduced. If the error estimate is acceptably small, the
663 * step is performed, success is returned and dt might be increased to make
664 * the steps as large as possible. This method also updates t if a step is
667 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
668 * Simple System concept.
669 * \param in The state of the ODE which should be solved.
670 * \param t The value of the time. Updated if the step is successful.
671 * \param out Used to store the result of the step.
672 * \param dt The step size. Updated.
673 * \return success if the step was accepted, fail otherwise.
675 template< class System , class StateIn , class StateOut >
676 typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type
677 try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
679 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 )
681 initialize( system , in , t );
683 return try_step( system , in , m_dxdt.m_v , t , out , dt );
688 * Version 3 : try_step( sys , x , dxdt , t , dt )
690 * This version does not solve the forwarding problem, boost::range can not be used.
693 * \brief Tries to perform one step.
695 * This method tries to do one step with step size dt. If the error estimate
696 * is to large, the step is rejected and the method returns fail and the
697 * step size dt is reduced. If the error estimate is acceptably small, the
698 * step is performed, success is returned and dt might be increased to make
699 * the steps as large as possible. This method also updates t if a step is
702 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
703 * Simple System concept.
704 * \param x The state of the ODE which should be solved. Overwritten if
705 * the step is successful.
706 * \param dxdt The derivative of state.
707 * \param t The value of the time. Updated if the step is successful.
708 * \param dt The step size. Updated.
709 * \return success if the step was accepted, fail otherwise.
711 template< class System , class StateInOut , class DerivInOut >
712 controlled_step_result try_step( System system , StateInOut &x , DerivInOut &dxdt , time_type &t , time_type &dt )
714 m_xnew_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_xnew_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
715 m_dxdt_new_resizer.adjust_size( x , detail::bind( &controlled_runge_kutta::template resize_m_dxdt_new_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
716 controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , m_dxdtnew.m_v , dt );
719 boost::numeric::odeint::copy( m_xnew.m_v , x );
720 boost::numeric::odeint::copy( m_dxdtnew.m_v , dxdt );
727 * Version 4 : try_step( sys , in , dxdt_in , t , out , dxdt_out , dt )
729 * This version does not solve the forwarding problem, boost::range can not be used.
732 * \brief Tries to perform one step.
734 * This method tries to do one step with step size dt. If the error estimate
735 * is to large, the step is rejected and the method returns fail and the
736 * step size dt is reduced. If the error estimate is acceptably small, the
737 * step is performed, success is returned and dt might be increased to make
738 * the steps as large as possible. This method also updates t if a step is
741 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
742 * Simple System concept.
743 * \param in The state of the ODE which should be solved.
744 * \param dxdt The derivative of state.
745 * \param t The value of the time. Updated if the step is successful.
746 * \param out Used to store the result of the step.
747 * \param dt The step size. Updated.
748 * \return success if the step was accepted, fail otherwise.
750 template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
751 controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t ,
752 StateOut &out , DerivOut &dxdt_out , time_type &dt )
754 if( !m_step_adjuster.check_step_size_limit(dt) )
756 // given dt was above step size limit - adjust and return fail;
757 dt = m_step_adjuster.get_max_dt();
761 m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
763 //fsal: m_stepper.get_dxdt( dxdt );
764 //fsal: m_stepper.do_step( sys , x , dxdt , t , dt , m_x_err );
765 m_stepper.do_step( system , in , dxdt_in , t , out , dxdt_out , dt , m_xerr.m_v );
767 // this potentially overwrites m_x_err! (standard_error_checker does, at least)
768 value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt_in , m_xerr.m_v , dt );
770 if( max_rel_err > 1.0 )
772 // error too big, decrease step size and reject this step
773 dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
776 // otherwise, increase step size and accept
778 dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
784 * \brief Resets the internal state of the underlying FSAL stepper.
792 * \brief Initializes the internal state storing an internal copy of the derivative.
794 * \param deriv The initial derivative of the ODE.
796 template< class DerivIn >
797 void initialize( const DerivIn &deriv )
799 boost::numeric::odeint::copy( deriv , m_dxdt.m_v );
800 m_first_call = false;
804 * \brief Initializes the internal state storing an internal copy of the derivative.
806 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the
807 * Simple System concept.
808 * \param x The initial state of the ODE which should be solved.
809 * \param t The initial time.
811 template< class System , class StateIn >
812 void initialize( System system , const StateIn &x , time_type t )
814 typename odeint::unwrap_reference< System >::type &sys = system;
815 sys( x , m_dxdt.m_v , t );
816 m_first_call = false;
820 * \brief Returns true if the stepper has been initialized, false otherwise.
822 * \return true, if the stepper has been initialized, false otherwise.
824 bool is_initialized( void ) const
826 return ! m_first_call;
831 * \brief Adjust the size of all temporaries in the stepper manually.
832 * \param x A state from which the size of the temporaries to be resized is deduced.
834 template< class StateType >
835 void adjust_size( const StateType &x )
837 resize_m_xerr_impl( x );
838 resize_m_dxdt_impl( x );
839 resize_m_dxdt_new_impl( x );
840 resize_m_xnew_impl( x );
845 * \brief Returns the instance of the underlying stepper.
846 * \returns The instance of the underlying stepper.
848 stepper_type& stepper( void )
854 * \brief Returns the instance of the underlying stepper.
855 * \returns The instance of the underlying stepper.
857 const stepper_type& stepper( void ) const
867 template< class StateIn >
868 bool resize_m_xerr_impl( const StateIn &x )
870 return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() );
873 template< class StateIn >
874 bool resize_m_dxdt_impl( const StateIn &x )
876 return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
879 template< class StateIn >
880 bool resize_m_dxdt_new_impl( const StateIn &x )
882 return adjust_size_by_resizeability( m_dxdtnew , x , typename is_resizeable<deriv_type>::type() );
885 template< class StateIn >
886 bool resize_m_xnew_impl( const StateIn &x )
888 return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
892 template< class System , class StateInOut >
893 controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
895 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 )
897 initialize( system , x , t );
899 return try_step( system , x , m_dxdt.m_v , t , dt );
903 stepper_type m_stepper;
904 error_checker_type m_error_checker;
905 step_adjuster_type m_step_adjuster;
907 resizer_type m_dxdt_resizer;
908 resizer_type m_xerr_resizer;
909 resizer_type m_xnew_resizer;
910 resizer_type m_dxdt_new_resizer;
912 wrapped_deriv_type m_dxdt;
913 wrapped_state_type m_xerr;
914 wrapped_state_type m_xnew;
915 wrapped_deriv_type m_dxdtnew;
920 /********** DOXYGEN **********/
922 /**** DEFAULT ERROR CHECKER ****/
925 * \class default_error_checker
926 * \brief The default error checker to be used with Runge-Kutta error steppers
928 * This class provides the default mechanism to compare the error estimates
929 * reported by Runge-Kutta error steppers with user defined error bounds.
930 * It is used by the controlled_runge_kutta steppers.
932 * \tparam Value The value type.
933 * \tparam Time The time type.
934 * \tparam Algebra The algebra type.
935 * \tparam Operations The operations type.
939 * \fn default_error_checker( value_type eps_abs , value_type eps_rel , value_type a_x , value_type a_dxdt ,
941 * \brief Constructs the error checker.
943 * The error is calculated as follows: ????
945 * \param eps_abs Absolute tolerance level.
946 * \param eps_rel Relative tolerance level.
947 * \param a_x Factor for the weight of the state.
948 * \param a_dxdt Factor for the weight of the derivative.
949 * \param max_dt Maximum allowed step size.
953 * \fn error( const State &x_old , const Deriv &dxdt_old , Err &x_err , time_type dt ) const
954 * \brief Calculates the error level.
956 * If the returned error level is greater than 1, the estimated error was
957 * larger than the permitted error bounds and the step should be repeated
958 * with a smaller step size.
960 * \param x_old State at the beginning of the step.
961 * \param dxdt_old Derivative at the beginning of the step.
962 * \param x_err Error estimate.
963 * \param dt Time step.
968 * \fn error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , time_type dt ) const
969 * \brief Calculates the error level using a given algebra.
971 * If the returned error level is greater than 1, the estimated error was
972 * larger than the permitted error bounds and the step should be repeated
973 * with a smaller step size.
975 * \param algebra The algebra used for calculation of the error.
976 * \param x_old State at the beginning of the step.
977 * \param dxdt_old Derivative at the beginning of the step.
978 * \param x_err Error estimate.
979 * \param dt Time step.
984 * \fn time_type decrease_step(const time_type dt, const value_type error, const int error_order)
985 * \brief Returns a decreased step size based on the given error and order
987 * Calculates a new smaller step size based on the given error and its order.
989 * \param dt The old step size.
990 * \param error The computed error estimate.
991 * \param error_order The error order of the stepper.
992 * \return dt_new The new, reduced step size.
996 * \fn time_type increase_step(const time_type dt, const value_type error, const int error_order)
997 * \brief Returns an increased step size based on the given error and order.
999 * Calculates a new bigger step size based on the given error and its order. If max_dt != 0, the
1000 * new step size is limited to max_dt.
1002 * \param dt The old step size.
1003 * \param error The computed error estimate.
1004 * \param error_order The order of the stepper.
1005 * \return dt_new The new, increased step size.
1014 #endif // BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_RUNGE_KUTTA_HPP_INCLUDED