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 rosenbrock4_controller(value_type atol=1.0e-6,value_type rtol=1.0e-6,const stepper_type & stepper=stepper_type ())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 rosenbrock4_controller(value_type atol,value_type rtol,time_type max_dt,const stepper_type & stepper=stepper_type ())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 error(const state_type & x,const state_type & xold,const state_type & xerr)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 last_error(void) const90 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 try_step(System sys,state_type & x,time_type & t,time_type & dt)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 try_step(System sys,const state_type & x,time_type & t,state_type & xout,time_type & dt)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 } else { 167 dt = dt_new; 168 } 169 m_last_rejected = false; 170 return success; 171 } 172 else 173 { 174 dt = dt_new; 175 m_last_rejected = true; 176 return fail; 177 } 178 } 179 180 181 template< class StateType > adjust_size(const StateType & x)182 void adjust_size( const StateType &x ) 183 { 184 resize_m_xerr( x ); 185 resize_m_xnew( x ); 186 } 187 188 189 stepper(void)190 stepper_type& stepper( void ) 191 { 192 return m_stepper; 193 } 194 stepper(void) const195 const stepper_type& stepper( void ) const 196 { 197 return m_stepper; 198 } 199 200 201 202 203 private: 204 205 template< class StateIn > resize_m_xerr(const StateIn & x)206 bool resize_m_xerr( const StateIn &x ) 207 { 208 return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() ); 209 } 210 211 template< class StateIn > resize_m_xnew(const StateIn & x)212 bool resize_m_xnew( const StateIn &x ) 213 { 214 return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() ); 215 } 216 217 218 stepper_type m_stepper; 219 resizer_type m_xerr_resizer; 220 resizer_type m_xnew_resizer; 221 wrapped_state_type m_xerr; 222 wrapped_state_type m_xnew; 223 value_type m_atol , m_rtol; 224 time_type m_max_dt; 225 bool m_first_step; 226 value_type m_err_old , m_dt_old; 227 bool m_last_rejected; 228 }; 229 230 231 232 233 234 235 } // namespace odeint 236 } // namespace numeric 237 } // namespace boost 238 239 240 #endif // BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED 241