ceph/src/boost/libs/numeric/odeint/include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp

   1 /*
   2  [auto_generated]
   3  boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
   4
   5  [begin_description]
   6  Controller for the Rosenbrock4 method.
   7  [end_description]
   8
   9  Copyright 2011-2012 Karsten Ahnert
  10  Copyright 2011-2012 Mario Mulansky
  11  Copyright 2012 Christoph Koke
  12
  13  Distributed under the Boost Software License, Version 1.0.
  14  (See accompanying file LICENSE_1_0.txt or
  15  copy at http://www.boost.org/LICENSE_1_0.txt)
  16  */
  17
  18
  19 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED
  20 #define BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED
  21
  22 #include <boost/config.hpp>
  23 #include <boost/numeric/odeint/util/bind.hpp>
  24
  25 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
  26 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
  27
  28 #include <boost/numeric/odeint/util/copy.hpp>
  29 #include <boost/numeric/odeint/util/is_resizeable.hpp>
  30 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
  31
  32 #include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
  33
  34 namespace boost {
  35 namespace numeric {
  36 namespace odeint {
  37
  38 template< class Stepper >
  39 class rosenbrock4_controller
  40 {
  41 private:
  42
  43
  44 public:
  45
  46     typedef Stepper stepper_type;
  47     typedef typename stepper_type::value_type value_type;
  48     typedef typename stepper_type::state_type state_type;
  49     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
  50     typedef typename stepper_type::time_type time_type;
  51     typedef typename stepper_type::deriv_type deriv_type;
  52     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
  53     typedef typename stepper_type::resizer_type resizer_type;
  54     typedef controlled_stepper_tag stepper_category;
  55
  56     typedef rosenbrock4_controller< Stepper > controller_type;
  57
  58
  59     rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 ,
  60                             const stepper_type &stepper = stepper_type() )
  61         : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) ,
  62           m_max_dt( static_cast<time_type>(0) ) ,
  63           m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
  64           m_last_rejected( false )
  65     { }
  66
  67     rosenbrock4_controller( value_type atol, value_type rtol, time_type max_dt,
  68                             const stepper_type &stepper = stepper_type() )
  69             : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) , m_max_dt( max_dt ) ,
  70               m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
  71               m_last_rejected( false )
  72     { }
  73
  74     value_type error( const state_type &x , const state_type &xold , const state_type &xerr )
  75     {
  76         BOOST_USING_STD_MAX();
  77         using std::abs;
  78         using std::sqrt;
  79
  80         const size_t n = x.size();
  81         value_type err = 0.0 , sk = 0.0;
  82         for( size_t i=0 ; i<n ; ++i )
  83         {
  84             sk = m_atol + m_rtol * max BOOST_PREVENT_MACRO_SUBSTITUTION ( abs( xold[i] ) , abs( x[i] ) );
  85             err += xerr[i] * xerr[i] / sk / sk;
  86         }
  87         return sqrt( err / value_type( n ) );
  88     }
  89
  90     value_type last_error( void ) const
  91     {
  92         return m_err_old;
  93     }
  94
  95
  96
  97
  98     template< class System >
  99     boost::numeric::odeint::controlled_step_result
 100     try_step( System sys , state_type &x , time_type &t , time_type &dt )
 101     {
 102         m_xnew_resizer.adjust_size( x , detail::bind( &controller_type::template resize_m_xnew< state_type > , detail::ref( *this ) , detail::_1 ) );
 103         boost::numeric::odeint::controlled_step_result res = try_step( sys , x , t , m_xnew.m_v , dt );
 104         if( res == success )
 105         {
 106             boost::numeric::odeint::copy( m_xnew.m_v , x );
 107         }
 108         return res;
 109     }
 110
 111
 112     template< class System >
 113     boost::numeric::odeint::controlled_step_result
 114     try_step( System sys , const state_type &x , time_type &t , state_type &xout , time_type &dt )
 115     {
 116         if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
 117         {
 118             // given step size is bigger then max_dt
 119             // set limit and return fail
 120             dt = m_max_dt;
 121             return fail;
 122         }
 123
 124         BOOST_USING_STD_MIN();
 125         BOOST_USING_STD_MAX();
 126         using std::pow;
 127
 128         static const value_type safe = 0.9 , fac1 = 5.0 , fac2 = 1.0 / 6.0;
 129
 130         m_xerr_resizer.adjust_size( x , detail::bind( &controller_type::template resize_m_xerr< state_type > , detail::ref( *this ) , detail::_1 ) );
 131
 132         m_stepper.do_step( sys , x , t , xout , dt , m_xerr.m_v );
 133         value_type err = error( xout , x , m_xerr.m_v );
 134
 135         value_type fac = max BOOST_PREVENT_MACRO_SUBSTITUTION (
 136             fac2 , min BOOST_PREVENT_MACRO_SUBSTITUTION (
 137                 fac1 ,
 138                 static_cast< value_type >( pow( err , 0.25 ) / safe ) ) );
 139         value_type dt_new = dt / fac;
 140         if ( err <= 1.0 )
 141         {
 142             if( m_first_step )
 143             {
 144                 m_first_step = false;
 145             }
 146             else
 147             {
 148                 value_type fac_pred = ( m_dt_old / dt ) * pow( err * err / m_err_old , 0.25 ) / safe;
 149                 fac_pred = max BOOST_PREVENT_MACRO_SUBSTITUTION (
 150                     fac2 , min BOOST_PREVENT_MACRO_SUBSTITUTION ( fac1 , fac_pred ) );
 151                 fac = max BOOST_PREVENT_MACRO_SUBSTITUTION ( fac , fac_pred );
 152                 dt_new = dt / fac;
 153             }
 154
 155             m_dt_old = dt;
 156             m_err_old = max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast< value_type >( 0.01 ) , err );
 157             if( m_last_rejected )
 158                 dt_new = ( dt >= 0.0 ?
 159                 min BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) :
 160                 max BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) );
 161             t += dt;
 162             // limit step size to max_dt
 163             if( m_max_dt != static_cast<time_type>(0) )
 164             {
 165                 dt = detail::min_abs(m_max_dt, dt_new);
 166             }
 167             m_last_rejected = false;
 168             return success;
 169         }
 170         else
 171         {
 172             dt = dt_new;
 173             m_last_rejected = true;
 174             return fail;
 175         }
 176     }
 177
 178
 179     template< class StateType >
 180     void adjust_size( const StateType &x )
 181     {
 182         resize_m_xerr( x );
 183         resize_m_xnew( x );
 184     }
 185
 186
 187
 188     stepper_type& stepper( void )
 189     {
 190         return m_stepper;
 191     }
 192
 193     const stepper_type& stepper( void ) const
 194     {
 195         return m_stepper;
 196     }
 197
 198
 199
 200
 201 private:
 202
 203     template< class StateIn >
 204     bool resize_m_xerr( const StateIn &x )
 205     {
 206         return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() );
 207     }
 208
 209     template< class StateIn >
 210     bool resize_m_xnew( const StateIn &x )
 211     {
 212         return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
 213     }
 214
 215
 216     stepper_type m_stepper;
 217     resizer_type m_xerr_resizer;
 218     resizer_type m_xnew_resizer;
 219     wrapped_state_type m_xerr;
 220     wrapped_state_type m_xnew;
 221     value_type m_atol , m_rtol;
 222     time_type m_max_dt;
 223     bool m_first_step;
 224     value_type m_err_old , m_dt_old;
 225     bool m_last_rejected;
 226 };
 227
 228
 229
 230
 231
 232
 233 } // namespace odeint
 234 } // namespace numeric
 235 } // namespace boost
 236
 237
 238 #endif // BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED