integrate/detail/integrate_const.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/integrate/detail/integrate_const.hpp
0004
0005  [begin_description]
0006  integrate const implementation
0007  [end_description]
0008
0009  Copyright 2012-2015 Mario Mulansky
0010  Copyright 2012 Christoph Koke
0011  Copyright 2012 Karsten Ahnert
0012
0013  Distributed under the Boost Software License, Version 1.0.
0014  (See accompanying file LICENSE_1_0.txt or
0015  copy at http://www.boost.org/LICENSE_1_0.txt)
0016  */
0017
0018 #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
0019 #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
0020
0021 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0022 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0023 #include <boost/numeric/odeint/util/unit_helper.hpp>
0024 #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
0025
0026 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
0027
0028 namespace boost {
0029 namespace numeric {
0030 namespace odeint {
0031 namespace detail {
0032
0033 // forward declaration
0034 template< class Stepper , class System , class State , class Time , class Observer >
0035 size_t integrate_adaptive(
0036         Stepper stepper , System system , State &start_state ,
0037         Time &start_time , Time end_time , Time &dt ,
0038         Observer observer , controlled_stepper_tag
0039 );
0040
0041
0042 template< class Stepper , class System , class State , class Time , class Observer >
0043 size_t integrate_const(
0044         Stepper stepper , System system , State &start_state ,
0045         Time start_time , Time end_time , Time dt ,
0046         Observer observer , stepper_tag
0047 )
0048 {
0049
0050     typename odeint::unwrap_reference< Observer >::type &obs = observer;
0051     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
0052
0053     Time time = start_time;
0054     int step = 0;
0055     // cast time+dt explicitely in case of expression templates (e.g. multiprecision)
0056     while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
0057     {
0058         obs( start_state , time );
0059         st.do_step( system , start_state , time , dt );
0060         // direct computation of the time avoids error propagation happening when using time += dt
0061         // we need clumsy type analysis to get boost units working here
0062         ++step;
0063         time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * dt;
0064     }
0065     obs( start_state , time );
0066
0067     return step;
0068 }
0069
0070
0071
0072 template< class Stepper , class System , class State , class Time , class Observer >
0073 size_t integrate_const(
0074         Stepper stepper , System system , State &start_state ,
0075         Time start_time , Time end_time , Time dt ,
0076         Observer observer , controlled_stepper_tag
0077 )
0078 {
0079     typename odeint::unwrap_reference< Observer >::type &obs = observer;
0080
0081     Time time = start_time;
0082     const Time time_step = dt;
0083     int real_steps = 0;
0084     int step = 0;
0085
0086     while( less_eq_with_sign( static_cast<Time>(time+time_step) , end_time , dt ) )
0087     {
0088         obs( start_state , time );
0089         // integrate_adaptive_checked uses the given checker to throw if an overflow occurs
0090         real_steps += detail::integrate_adaptive(stepper, system, start_state, time,
0091                                                  static_cast<Time>(time + time_step), dt,
0092                                                  null_observer(), controlled_stepper_tag());
0093         // direct computation of the time avoids error propagation happening when using time += dt
0094         // we need clumsy type analysis to get boost units working here
0095         step++;
0096         time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * time_step;
0097     }
0098     obs( start_state , time );
0099
0100     return real_steps;
0101 }
0102
0103
0104 template< class Stepper , class System , class State , class Time , class Observer >
0105 size_t integrate_const(
0106         Stepper stepper , System system , State &start_state ,
0107         Time start_time , Time end_time , Time dt ,
0108         Observer observer , dense_output_stepper_tag
0109 )
0110 {
0111     typename odeint::unwrap_reference< Observer >::type &obs = observer;
0112     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
0113
0114     Time time = start_time;
0115
0116     st.initialize( start_state , time , dt );
0117     obs( start_state , time );
0118     time += dt;
0119
0120     int obs_step( 1 );
0121     int real_step( 0 );
0122
0123     while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
0124     {
0125         while( less_eq_with_sign( time , st.current_time() , dt ) )
0126         {
0127             st.calc_state( time , start_state );
0128             obs( start_state , time );
0129             ++obs_step;
0130             // direct computation of the time avoids error propagation happening when using time += dt
0131             // we need clumsy type analysis to get boost units working here
0132             time = start_time + static_cast< typename unit_value_type<Time>::type >(obs_step) * dt;
0133         }
0134         // we have not reached the end, do another real step
0135         if( less_with_sign( static_cast<Time>(st.current_time()+st.current_time_step()) ,
0136                             end_time ,
0137                             st.current_time_step() ) )
0138         {
0139             while( less_eq_with_sign( st.current_time() , time , dt ) )
0140             {
0141                 st.do_step( system );
0142                 ++real_step;
0143             }
0144         }
0145         else if( less_with_sign( st.current_time() , end_time , st.current_time_step() ) )
0146         { // do the last step ending exactly on the end point
0147             st.initialize( st.current_state() , st.current_time() , end_time - st.current_time() );
0148             st.do_step( system );
0149             ++real_step;
0150         }
0151
0152     }
0153     // last observation, if we are still in observation interval
0154     // might happen due to finite precision problems when computing the the time
0155     if( less_eq_with_sign( time , end_time , dt ) )
0156     {
0157         st.calc_state( time , start_state );
0158         obs( start_state , time );
0159     }
0160
0161     return real_step;
0162 }
0163
0164
0165 } } } }
0166
0167 #endif