• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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