Chris@16: /*
Chris@16:  [auto_generated]
Chris@16:  boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
Chris@16: 
Chris@16:  [begin_description]
Chris@16:  Controller for the Rosenbrock4 method.
Chris@16:  [end_description]
Chris@16: 
Chris@101:  Copyright 2011-2012 Karsten Ahnert
Chris@101:  Copyright 2011-2012 Mario Mulansky
Chris@101:  Copyright 2012 Christoph Koke
Chris@16: 
Chris@16:  Distributed under the Boost Software License, Version 1.0.
Chris@16:  (See accompanying file LICENSE_1_0.txt or
Chris@16:  copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16:  */
Chris@16: 
Chris@16: 
Chris@16: #ifndef BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED
Chris@16: #define BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED
Chris@16: 
Chris@16: #include <boost/config.hpp>
Chris@16: #include <boost/numeric/odeint/util/bind.hpp>
Chris@16: 
Chris@16: #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
Chris@16: #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
Chris@16: 
Chris@16: #include <boost/numeric/odeint/util/copy.hpp>
Chris@16: #include <boost/numeric/odeint/util/is_resizeable.hpp>
Chris@16: 
Chris@16: #include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
Chris@16: 
Chris@16: namespace boost {
Chris@16: namespace numeric {
Chris@16: namespace odeint {
Chris@16: 
Chris@16: template< class Stepper >
Chris@16: class rosenbrock4_controller
Chris@16: {
Chris@16: private:
Chris@16: 
Chris@16: 
Chris@16: public:
Chris@16: 
Chris@16:     typedef Stepper stepper_type;
Chris@16:     typedef typename stepper_type::value_type value_type;
Chris@16:     typedef typename stepper_type::state_type state_type;
Chris@16:     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
Chris@16:     typedef typename stepper_type::time_type time_type;
Chris@16:     typedef typename stepper_type::deriv_type deriv_type;
Chris@16:     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
Chris@16:     typedef typename stepper_type::resizer_type resizer_type;
Chris@16:     typedef controlled_stepper_tag stepper_category;
Chris@16: 
Chris@16:     typedef rosenbrock4_controller< Stepper > controller_type;
Chris@16: 
Chris@16: 
Chris@16:     rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 , const stepper_type &stepper = stepper_type() )
Chris@16:     : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) ,
Chris@16:       m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
Chris@16:       m_last_rejected( false )
Chris@16:     { }
Chris@16: 
Chris@16: 
Chris@16:     value_type error( const state_type &x , const state_type &xold , const state_type &xerr )
Chris@16:     {
Chris@16:         BOOST_USING_STD_MAX();
Chris@16:         using std::abs;
Chris@101:         using std::sqrt;
Chris@16:         
Chris@16:         const size_t n = x.size();
Chris@16:         value_type err = 0.0 , sk = 0.0;
Chris@16:         for( size_t i=0 ; i<n ; ++i )
Chris@16:         {
Chris@16:             sk = m_atol + m_rtol * max BOOST_PREVENT_MACRO_SUBSTITUTION ( abs( xold[i] ) , abs( x[i] ) );
Chris@16:             err += xerr[i] * xerr[i] / sk / sk;
Chris@16:         }
Chris@16:         return sqrt( err / value_type( n ) );
Chris@16:     }
Chris@16: 
Chris@16:     value_type last_error( void ) const
Chris@16:     {
Chris@16:         return m_err_old;
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16:     template< class System >
Chris@16:     boost::numeric::odeint::controlled_step_result
Chris@16:     try_step( System sys , state_type &x , time_type &t , time_type &dt )
Chris@16:     {
Chris@16:         m_xnew_resizer.adjust_size( x , detail::bind( &controller_type::template resize_m_xnew< state_type > , detail::ref( *this ) , detail::_1 ) );
Chris@16:         boost::numeric::odeint::controlled_step_result res = try_step( sys , x , t , m_xnew.m_v , dt );
Chris@16:         if( res == success )
Chris@16:         {
Chris@16:             boost::numeric::odeint::copy( m_xnew.m_v , x );
Chris@16:         }
Chris@16:         return res;
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16:     template< class System >
Chris@16:     boost::numeric::odeint::controlled_step_result
Chris@16:     try_step( System sys , const state_type &x , time_type &t , state_type &xout , time_type &dt )
Chris@16:     {
Chris@16:         BOOST_USING_STD_MIN();
Chris@16:         BOOST_USING_STD_MAX();
Chris@16:         using std::pow;
Chris@16: 
Chris@16:         static const value_type safe = 0.9 , fac1 = 5.0 , fac2 = 1.0 / 6.0;
Chris@16: 
Chris@16:         m_xerr_resizer.adjust_size( x , detail::bind( &controller_type::template resize_m_xerr< state_type > , detail::ref( *this ) , detail::_1 ) );
Chris@16: 
Chris@16:         m_stepper.do_step( sys , x , t , xout , dt , m_xerr.m_v );
Chris@16:         value_type err = error( xout , x , m_xerr.m_v );
Chris@16: 
Chris@101:         value_type fac = max BOOST_PREVENT_MACRO_SUBSTITUTION (
Chris@101:             fac2 , min BOOST_PREVENT_MACRO_SUBSTITUTION (
Chris@101:                 fac1 ,
Chris@101:                 static_cast< value_type >( pow( err , 0.25 ) / safe ) ) );
Chris@16:         value_type dt_new = dt / fac;
Chris@16:         if ( err <= 1.0 )
Chris@16:         {
Chris@16:             if( m_first_step )
Chris@16:             {
Chris@16:                 m_first_step = false;
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                 value_type fac_pred = ( m_dt_old / dt ) * pow( err * err / m_err_old , 0.25 ) / safe;
Chris@101:                 fac_pred = max BOOST_PREVENT_MACRO_SUBSTITUTION (
Chris@101:                     fac2 , min BOOST_PREVENT_MACRO_SUBSTITUTION ( fac1 , fac_pred ) );
Chris@16:                 fac = max BOOST_PREVENT_MACRO_SUBSTITUTION ( fac , fac_pred );
Chris@16:                 dt_new = dt / fac;
Chris@16:             }
Chris@16: 
Chris@16:             m_dt_old = dt;
Chris@101:             m_err_old = max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast< value_type >( 0.01 ) , err );
Chris@16:             if( m_last_rejected )
Chris@101:                 dt_new = ( dt >= 0.0 ?
Chris@101:                 min BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) :
Chris@101:                 max BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) );
Chris@16:             t += dt;
Chris@16:             dt = dt_new;
Chris@16:             m_last_rejected = false;
Chris@16:             return success;
Chris@16:         }
Chris@16:         else
Chris@16:         {
Chris@16:             dt = dt_new;
Chris@16:             m_last_rejected = true;
Chris@16:             return fail;
Chris@16:         }
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16:     template< class StateType >
Chris@16:     void adjust_size( const StateType &x )
Chris@16:     {
Chris@16:         resize_m_xerr( x );
Chris@16:         resize_m_xnew( x );
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16:     stepper_type& stepper( void )
Chris@16:     {
Chris@16:         return m_stepper;
Chris@16:     }
Chris@16: 
Chris@16:     const stepper_type& stepper( void ) const
Chris@16:     {
Chris@16:         return m_stepper;
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: private:
Chris@16: 
Chris@16:     template< class StateIn >
Chris@16:     bool resize_m_xerr( const StateIn &x )
Chris@16:     {
Chris@16:         return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() );
Chris@16:     }
Chris@16: 
Chris@16:     template< class StateIn >
Chris@16:     bool resize_m_xnew( const StateIn &x )
Chris@16:     {
Chris@16:         return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16:     stepper_type m_stepper;
Chris@16:     resizer_type m_xerr_resizer;
Chris@16:     resizer_type m_xnew_resizer;
Chris@16:     wrapped_state_type m_xerr;
Chris@16:     wrapped_state_type m_xnew;
Chris@16:     value_type m_atol , m_rtol;
Chris@16:     bool m_first_step;
Chris@16:     value_type m_err_old , m_dt_old;
Chris@16:     bool m_last_rejected;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: } // namespace odeint
Chris@16: } // namespace numeric
Chris@16: } // namespace boost
Chris@16: 
Chris@16: 
Chris@16: #endif // BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED