1 /* 2 [auto_generated] 3 boost/numeric/odeint/stepper/adams_bashforth.hpp 4 5 [begin_description] 6 Implementaton of the Adam-Bashforth method a multistep method used for the predictor step in the 7 Adams-Bashforth-Moulton method. 8 [end_description] 9 10 Copyright 2011-2013 Karsten Ahnert 11 Copyright 2011-2013 Mario Mulansky 12 Copyright 2012 Christoph Koke 13 Copyright 2013 Pascal Germroth 14 15 Distributed under the Boost Software License, Version 1.0. 16 (See accompanying file LICENSE_1_0.txt or 17 copy at http://www.boost.org/LICENSE_1_0.txt) 18 */ 19 20 21 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_ADAMS_BASHFORTH_HPP_INCLUDED 22 #define BOOST_NUMERIC_ODEINT_STEPPER_ADAMS_BASHFORTH_HPP_INCLUDED 23 24 #include <boost/static_assert.hpp> 25 26 #include <boost/numeric/odeint/util/bind.hpp> 27 #include <boost/numeric/odeint/util/unwrap_reference.hpp> 28 29 #include <boost/numeric/odeint/algebra/range_algebra.hpp> 30 #include <boost/numeric/odeint/algebra/default_operations.hpp> 31 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> 32 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> 33 34 #include <boost/numeric/odeint/util/state_wrapper.hpp> 35 #include <boost/numeric/odeint/util/is_resizeable.hpp> 36 #include <boost/numeric/odeint/util/resizer.hpp> 37 38 #include <boost/numeric/odeint/stepper/stepper_categories.hpp> 39 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp> 40 #include <boost/numeric/odeint/stepper/extrapolation_stepper.hpp> 41 42 #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp> 43 44 #include <boost/numeric/odeint/stepper/detail/adams_bashforth_coefficients.hpp> 45 #include <boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp> 46 #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp> 47 48 #include <boost/mpl/arithmetic.hpp> 49 #include <boost/mpl/min_max.hpp> 50 #include <boost/mpl/equal_to.hpp> 51 52 namespace mpl = boost::mpl; 53 54 55 namespace boost { 56 namespace numeric { 57 namespace odeint { 58 59 using mpl::int_; 60 61 /* if N >= 4, returns the smallest even number > N, otherwise returns 4 */ 62 template < int N > 63 struct order_helper 64 : mpl::max< typename mpl::eval_if< 65 mpl::equal_to< mpl::modulus< int_< N >, int_< 2 > >, 66 int_< 0 > >, 67 int_< N >, int_< N + 1 > >::type, 68 int_< 4 > >::type 69 { }; 70 71 template< 72 size_t Steps , 73 class State , 74 class Value = double , 75 class Deriv = State , 76 class Time = Value , 77 class Algebra = typename algebra_dispatcher< State >::algebra_type , 78 class Operations = typename operations_dispatcher< State >::operations_type , 79 class Resizer = initially_resizer , 80 class InitializingStepper = extrapolation_stepper< order_helper<Steps>::value, 81 State, Value, Deriv, Time, 82 Algebra, Operations, Resizer > 83 > 84 class adams_bashforth : public algebra_stepper_base< Algebra , Operations > 85 { 86 87 #ifndef DOXYGEN_SKIP 88 BOOST_STATIC_ASSERT(( Steps > 0 )); 89 BOOST_STATIC_ASSERT(( Steps < 9 )); 90 #endif 91 92 public : 93 94 typedef State state_type; 95 typedef state_wrapper< state_type > wrapped_state_type; 96 typedef Value value_type; 97 typedef Deriv deriv_type; 98 typedef state_wrapper< deriv_type > wrapped_deriv_type; 99 typedef Time time_type; 100 typedef Resizer resizer_type; 101 typedef stepper_tag stepper_category; 102 103 typedef InitializingStepper initializing_stepper_type; 104 105 typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type; 106 typedef typename algebra_stepper_base_type::algebra_type algebra_type; 107 typedef typename algebra_stepper_base_type::operations_type operations_type; 108 #ifndef DOXYGEN_SKIP 109 typedef adams_bashforth< Steps , State , Value , Deriv , Time , Algebra , Operations , Resizer , InitializingStepper > stepper_type; 110 #endif 111 static const size_t steps = Steps; 112 113 114 115 typedef unsigned short order_type; 116 static const order_type order_value = steps; 117 118 typedef detail::rotating_buffer< wrapped_deriv_type , steps > step_storage_type; 119 120 121 order(void) const122 order_type order( void ) const { return order_value; } 123 adams_bashforth(const algebra_type & algebra=algebra_type ())124 adams_bashforth( const algebra_type &algebra = algebra_type() ) 125 : algebra_stepper_base_type( algebra ) , 126 m_step_storage() , m_resizer() , m_coefficients() , 127 m_steps_initialized( 0 ) , m_initializing_stepper() 128 { } 129 130 131 132 /* 133 * Version 1 : do_step( system , x , t , dt ); 134 * 135 * solves the forwarding problem 136 */ 137 template< class System , class StateInOut > do_step(System system,StateInOut & x,time_type t,time_type dt)138 void do_step( System system , StateInOut &x , time_type t , time_type dt ) 139 { 140 do_step( system , x , t , x , dt ); 141 } 142 143 /** 144 * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateInOut. 145 */ 146 template< class System , class StateInOut > do_step(System system,const StateInOut & x,time_type t,time_type dt)147 void do_step( System system , const StateInOut &x , time_type t , time_type dt ) 148 { 149 do_step( system , x , t , x , dt ); 150 } 151 152 153 154 /* 155 * Version 2 : do_step( system , in , t , out , dt ); 156 * 157 * solves the forwarding problem 158 */ 159 160 template< class System , class StateIn , class StateOut > do_step(System system,const StateIn & in,time_type t,StateOut & out,time_type dt)161 void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt ) 162 { 163 do_step_impl( system , in , t , out , dt ); 164 } 165 166 /** 167 * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateOut. 168 */ 169 template< class System , class StateIn , class StateOut > do_step(System system,const StateIn & in,time_type t,const StateOut & out,time_type dt)170 void do_step( System system , const StateIn &in , time_type t , const StateOut &out , time_type dt ) 171 { 172 do_step_impl( system , in , t , out , dt ); 173 } 174 175 176 template< class StateType > adjust_size(const StateType & x)177 void adjust_size( const StateType &x ) 178 { 179 resize_impl( x ); 180 } 181 step_storage(void) const182 const step_storage_type& step_storage( void ) const 183 { 184 return m_step_storage; 185 } 186 step_storage(void)187 step_storage_type& step_storage( void ) 188 { 189 return m_step_storage; 190 } 191 192 template< class ExplicitStepper , class System , class StateIn > initialize(ExplicitStepper explicit_stepper,System system,StateIn & x,time_type & t,time_type dt)193 void initialize( ExplicitStepper explicit_stepper , System system , StateIn &x , time_type &t , time_type dt ) 194 { 195 typename odeint::unwrap_reference< ExplicitStepper >::type &stepper = explicit_stepper; 196 typename odeint::unwrap_reference< System >::type &sys = system; 197 198 m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) ); 199 200 for( size_t i=0 ; i+1<steps ; ++i ) 201 { 202 if( i != 0 ) m_step_storage.rotate(); 203 sys( x , m_step_storage[0].m_v , t ); 204 stepper.do_step_dxdt_impl( system, x, m_step_storage[0].m_v, t, 205 dt ); 206 t += dt; 207 } 208 m_steps_initialized = steps; 209 } 210 211 template< class System , class StateIn > initialize(System system,StateIn & x,time_type & t,time_type dt)212 void initialize( System system , StateIn &x , time_type &t , time_type dt ) 213 { 214 initialize( detail::ref( m_initializing_stepper ) , system , x , t , dt ); 215 } 216 reset(void)217 void reset( void ) 218 { 219 m_steps_initialized = 0; 220 } 221 is_initialized(void) const222 bool is_initialized( void ) const 223 { 224 return m_steps_initialized >= ( steps - 1 ); 225 } 226 initializing_stepper(void) const227 const initializing_stepper_type& initializing_stepper( void ) const { return m_initializing_stepper; } 228 initializing_stepper(void)229 initializing_stepper_type& initializing_stepper( void ) { return m_initializing_stepper; } 230 231 private: 232 233 template< class System , class StateIn , class StateOut > do_step_impl(System system,const StateIn & in,time_type t,StateOut & out,time_type dt)234 void do_step_impl( System system , const StateIn &in , time_type t , StateOut &out , time_type dt ) 235 { 236 typename odeint::unwrap_reference< System >::type &sys = system; 237 if( m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) ) ) 238 { 239 m_steps_initialized = 0; 240 } 241 242 if( m_steps_initialized + 1 < steps ) 243 { 244 if( m_steps_initialized != 0 ) m_step_storage.rotate(); 245 sys( in , m_step_storage[0].m_v , t ); 246 m_initializing_stepper.do_step_dxdt_impl( 247 system, in, m_step_storage[0].m_v, t, out, dt ); 248 ++m_steps_initialized; 249 } 250 else 251 { 252 m_step_storage.rotate(); 253 sys( in , m_step_storage[0].m_v , t ); 254 detail::adams_bashforth_call_algebra< steps , algebra_type , operations_type >()( this->m_algebra , in , out , m_step_storage , m_coefficients , dt ); 255 } 256 } 257 258 259 template< class StateIn > resize_impl(const StateIn & x)260 bool resize_impl( const StateIn &x ) 261 { 262 bool resized( false ); 263 for( size_t i=0 ; i<steps ; ++i ) 264 { 265 resized |= adjust_size_by_resizeability( m_step_storage[i] , x , typename is_resizeable<deriv_type>::type() ); 266 } 267 return resized; 268 } 269 270 step_storage_type m_step_storage; 271 resizer_type m_resizer; 272 detail::adams_bashforth_coefficients< value_type , steps > m_coefficients; 273 size_t m_steps_initialized; 274 initializing_stepper_type m_initializing_stepper; 275 276 }; 277 278 279 /***** DOXYGEN *****/ 280 281 /** 282 * \class adams_bashforth 283 * \brief The Adams-Bashforth multistep algorithm. 284 * 285 * The Adams-Bashforth method is a multi-step algorithm with configurable step 286 * number. The step number is specified as template parameter Steps and it 287 * then uses the result from the previous Steps steps. See also 288 * <a href="http://en.wikipedia.org/wiki/Linear_multistep_method">en.wikipedia.org/wiki/Linear_multistep_method</a>. 289 * Currently, a maximum of Steps=8 is supported. 290 * The method is explicit and fulfills the Stepper concept. Step size control 291 * or continuous output are not provided. 292 * 293 * This class derives from algebra_base and inherits its interface via 294 * CRTP (current recurring template pattern). For more details see 295 * algebra_stepper_base. 296 * 297 * \tparam Steps The number of steps (maximal 8). 298 * \tparam State The state type. 299 * \tparam Value The value type. 300 * \tparam Deriv The type representing the time derivative of the state. 301 * \tparam Time The time representing the independent variable - the time. 302 * \tparam Algebra The algebra type. 303 * \tparam Operations The operations type. 304 * \tparam Resizer The resizer policy type. 305 * \tparam InitializingStepper The stepper for the first two steps. 306 */ 307 308 /** 309 * \fn adams_bashforth::adams_bashforth( const algebra_type &algebra ) 310 * \brief Constructs the adams_bashforth class. This constructor can be used as a default 311 * constructor if the algebra has a default constructor. 312 * \param algebra A copy of algebra is made and stored. 313 */ 314 315 /** 316 * \fn order_type adams_bashforth::order( void ) const 317 * \brief Returns the order of the algorithm, which is equal to the number of steps. 318 * \return order of the method. 319 */ 320 321 /** 322 * \fn void adams_bashforth::do_step( System system , StateInOut &x , time_type t , time_type dt ) 323 * \brief This method performs one step. It transforms the result in-place. 324 * 325 * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the 326 * Simple System concept. 327 * \param x The state of the ODE which should be solved. After calling do_step the result is updated in x. 328 * \param t The value of the time, at which the step should be performed. 329 * \param dt The step size. 330 */ 331 332 /** 333 * \fn void adams_bashforth::do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt ) 334 * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place. 335 * 336 * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the 337 * Simple System concept. 338 * \param in The state of the ODE which should be solved. in is not modified in this method 339 * \param t The value of the time, at which the step should be performed. 340 * \param out The result of the step is written in out. 341 * \param dt The step size. 342 */ 343 344 /** 345 * \fn void adams_bashforth::adjust_size( const StateType &x ) 346 * \brief Adjust the size of all temporaries in the stepper manually. 347 * \param x A state from which the size of the temporaries to be resized is deduced. 348 */ 349 350 351 /** 352 * \fn const step_storage_type& adams_bashforth::step_storage( void ) const 353 * \brief Returns the storage of intermediate results. 354 * \return The storage of intermediate results. 355 */ 356 357 /** 358 * \fn step_storage_type& adams_bashforth::step_storage( void ) 359 * \brief Returns the storage of intermediate results. 360 * \return The storage of intermediate results. 361 */ 362 363 /** 364 * \fn void adams_bashforth::initialize( ExplicitStepper explicit_stepper , System system , StateIn &x , time_type &t , time_type dt ) 365 * \brief Initialized the stepper. Does Steps-1 steps with the explicit_stepper to fill the buffer. 366 * \param explicit_stepper the stepper used to fill the buffer of previous step results 367 * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the 368 * Simple System concept. 369 * \param x The state of the ODE which should be solved. After calling do_step the result is updated in x. 370 * \param t The value of the time, at which the step should be performed. 371 * \param dt The step size. 372 */ 373 374 /** 375 * \fn void adams_bashforth::initialize( System system , StateIn &x , time_type &t , time_type dt ) 376 * \brief Initialized the stepper. Does Steps-1 steps with an internal instance of InitializingStepper to fill the buffer. 377 * \note The state x and time t are updated to the values after Steps-1 initial steps. 378 * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the 379 * Simple System concept. 380 * \param x The initial state of the ODE which should be solved, updated in this method. 381 * \param t The initial value of the time, updated in this method. 382 * \param dt The step size. 383 */ 384 385 /** 386 * \fn void adams_bashforth::reset( void ) 387 * \brief Resets the internal buffer of the stepper. 388 */ 389 390 /** 391 * \fn bool adams_bashforth::is_initialized( void ) const 392 * \brief Returns true if the stepper has been initialized. 393 * \return bool true if stepper is initialized, false otherwise 394 */ 395 396 /** 397 * \fn const initializing_stepper_type& adams_bashforth::initializing_stepper( void ) const 398 * \brief Returns the internal initializing stepper instance. 399 * \return initializing_stepper 400 */ 401 402 /** 403 * \fn const initializing_stepper_type& adams_bashforth::initializing_stepper( void ) const 404 * \brief Returns the internal initializing stepper instance. 405 * \return initializing_stepper 406 */ 407 408 /** 409 * \fn initializing_stepper_type& adams_bashforth::initializing_stepper( void ) 410 * \brief Returns the internal initializing stepper instance. 411 * \return initializing_stepper 412 */ 413 414 } // odeint 415 } // numeric 416 } // boost 417 418 419 420 #endif // BOOST_NUMERIC_ODEINT_STEPPER_ADAMS_BASHFORTH_HPP_INCLUDED 421