annotate DEPENDENCIES/generic/include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp @ 46:d572322e2efe

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