odeint/stepper/rosenbrock4_controller.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
0004
0005  [begin_description]
0006  Controller for the Rosenbrock4 method.
0007  [end_description]
0008
0009  Copyright 2011-2012 Karsten Ahnert
0010  Copyright 2011-2012 Mario Mulansky
0011  Copyright 2012 Christoph Koke
0012
0013  Distributed under the Boost Software License, Version 1.0.
0014  (See accompanying file LICENSE_1_0.txt or
0015  copy at http://www.boost.org/LICENSE_1_0.txt)
0016  */
0017
0018
0019 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED
0020 #define BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED
0021
0022 #include <boost/config.hpp>
0023 #include <boost/numeric/odeint/util/bind.hpp>
0024
0025 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
0026 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0027
0028 #include <boost/numeric/odeint/util/copy.hpp>
0029 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0030 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
0031
0032 #include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
0033
0034 namespace boost {
0035 namespace numeric {
0036 namespace odeint {
0037
0038 template< class Stepper >
0039 class rosenbrock4_controller
0040 {
0041 private:
0042
0043
0044 public:
0045
0046     typedef Stepper stepper_type;
0047     typedef typename stepper_type::value_type value_type;
0048     typedef typename stepper_type::state_type state_type;
0049     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
0050     typedef typename stepper_type::time_type time_type;
0051     typedef typename stepper_type::deriv_type deriv_type;
0052     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
0053     typedef typename stepper_type::resizer_type resizer_type;
0054     typedef controlled_stepper_tag stepper_category;
0055
0056     typedef rosenbrock4_controller< Stepper > controller_type;
0057
0058
0059     rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 ,
0060                             const stepper_type &stepper = stepper_type() )
0061         : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) ,
0062           m_max_dt( static_cast<time_type>(0) ) ,
0063           m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
0064           m_last_rejected( false )
0065     { }
0066
0067     rosenbrock4_controller( value_type atol, value_type rtol, time_type max_dt,
0068                             const stepper_type &stepper = stepper_type() )
0069             : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) , m_max_dt( max_dt ) ,
0070               m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
0071               m_last_rejected( false )
0072     { }
0073
0074     value_type error( const state_type &x , const state_type &xold , const state_type &xerr )
0075     {
0076         BOOST_USING_STD_MAX();
0077         using std::abs;
0078         using std::sqrt;
0079
0080         const size_t n = x.size();
0081         value_type err = 0.0 , sk = 0.0;
0082         for( size_t i=0 ; i<n ; ++i )
0083         {
0084             sk = m_atol + m_rtol * max BOOST_PREVENT_MACRO_SUBSTITUTION ( abs( xold[i] ) , abs( x[i] ) );
0085             err += xerr[i] * xerr[i] / sk / sk;
0086         }
0087         return sqrt( err / value_type( n ) );
0088     }
0089
0090     value_type last_error( void ) const
0091     {
0092         return m_err_old;
0093     }
0094
0095
0096
0097
0098     template< class System >
0099     boost::numeric::odeint::controlled_step_result
0100     try_step( System sys , state_type &x , time_type &t , time_type &dt )
0101     {
0102         m_xnew_resizer.adjust_size( x , detail::bind( &controller_type::template resize_m_xnew< state_type > , detail::ref( *this ) , detail::_1 ) );
0103         boost::numeric::odeint::controlled_step_result res = try_step( sys , x , t , m_xnew.m_v , dt );
0104         if( res == success )
0105         {
0106             boost::numeric::odeint::copy( m_xnew.m_v , x );
0107         }
0108         return res;
0109     }
0110
0111
0112     template< class System >
0113     boost::numeric::odeint::controlled_step_result
0114     try_step( System sys , const state_type &x , time_type &t , state_type &xout , time_type &dt )
0115     {
0116         if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
0117         {
0118             // given step size is bigger then max_dt
0119             // set limit and return fail
0120             dt = m_max_dt;
0121             return fail;
0122         }
0123
0124         BOOST_USING_STD_MIN();
0125         BOOST_USING_STD_MAX();
0126         using std::pow;
0127
0128         static const value_type safe = 0.9 , fac1 = 5.0 , fac2 = 1.0 / 6.0;
0129
0130         m_xerr_resizer.adjust_size( x , detail::bind( &controller_type::template resize_m_xerr< state_type > , detail::ref( *this ) , detail::_1 ) );
0131
0132         m_stepper.do_step( sys , x , t , xout , dt , m_xerr.m_v );
0133         value_type err = error( xout , x , m_xerr.m_v );
0134
0135         value_type fac = max BOOST_PREVENT_MACRO_SUBSTITUTION (
0136             fac2 , min BOOST_PREVENT_MACRO_SUBSTITUTION (
0137                 fac1 ,
0138                 static_cast< value_type >( pow( err , 0.25 ) / safe ) ) );
0139         value_type dt_new = dt / fac;
0140         if ( err <= 1.0 )
0141         {
0142             if( m_first_step )
0143             {
0144                 m_first_step = false;
0145             }
0146             else
0147             {
0148                 value_type fac_pred = ( m_dt_old / dt ) * pow( err * err / m_err_old , 0.25 ) / safe;
0149                 fac_pred = max BOOST_PREVENT_MACRO_SUBSTITUTION (
0150                     fac2 , min BOOST_PREVENT_MACRO_SUBSTITUTION ( fac1 , fac_pred ) );
0151                 fac = max BOOST_PREVENT_MACRO_SUBSTITUTION ( fac , fac_pred );
0152                 dt_new = dt / fac;
0153             }
0154
0155             m_dt_old = dt;
0156             m_err_old = max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast< value_type >( 0.01 ) , err );
0157             if( m_last_rejected )
0158                 dt_new = ( dt >= 0.0 ?
0159                 min BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) :
0160                 max BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) );
0161             t += dt;
0162             // limit step size to max_dt
0163             if( m_max_dt != static_cast<time_type>(0) )
0164             {
0165                 dt = detail::min_abs(m_max_dt, dt_new);
0166             } else {
0167                 dt = dt_new;
0168             }
0169             m_last_rejected = false;
0170             return success;
0171         }
0172         else
0173         {
0174             dt = dt_new;
0175             m_last_rejected = true;
0176             return fail;
0177         }
0178     }
0179
0180
0181     template< class StateType >
0182     void adjust_size( const StateType &x )
0183     {
0184         resize_m_xerr( x );
0185         resize_m_xnew( x );
0186     }
0187
0188
0189
0190     stepper_type& stepper( void )
0191     {
0192         return m_stepper;
0193     }
0194
0195     const stepper_type& stepper( void ) const
0196     {
0197         return m_stepper;
0198     }
0199
0200
0201
0202
0203 protected:
0204
0205     template< class StateIn >
0206     bool resize_m_xerr( const StateIn &x )
0207     {
0208         return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() );
0209     }
0210
0211     template< class StateIn >
0212     bool resize_m_xnew( const StateIn &x )
0213     {
0214         return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
0215     }
0216
0217
0218     stepper_type m_stepper;
0219     resizer_type m_xerr_resizer;
0220     resizer_type m_xnew_resizer;
0221     wrapped_state_type m_xerr;
0222     wrapped_state_type m_xnew;
0223     value_type m_atol , m_rtol;
0224     time_type m_max_dt;
0225     bool m_first_step;
0226     value_type m_err_old , m_dt_old;
0227     bool m_last_rejected;
0228 };
0229
0230
0231
0232
0233
0234
0235 } // namespace odeint
0236 } // namespace numeric
0237 } // namespace boost
0238
0239
0240 #endif // BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED