integrate/detail/integrate_times.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/integrate/detail/integrate_times.hpp
0004
0005  [begin_description]
0006  Default integrate times implementation.
0007  [end_description]
0008
0009  Copyright 2011-2015 Mario Mulansky
0010  Copyright 2012 Karsten Ahnert
0011  Copyright 2012 Christoph Koke
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
0019 #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
0020 #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
0021
0022 #include <stdexcept>
0023
0024 #include <boost/config.hpp>
0025 #include <boost/throw_exception.hpp>
0026 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0027 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
0028 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
0029 #include <boost/numeric/odeint/integrate/max_step_checker.hpp>
0030
0031
0032 namespace boost {
0033 namespace numeric {
0034 namespace odeint {
0035 namespace detail {
0036
0037
0038
0039 /*
0040  * integrate_times for simple stepper
0041  */
0042 template<class Stepper, class System, class State, class TimeIterator, class Time, class Observer>
0043 size_t integrate_times(
0044         Stepper stepper , System system , State &start_state ,
0045         TimeIterator start_time , TimeIterator end_time , Time dt ,
0046         Observer observer , stepper_tag
0047 )
0048 {
0049     typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
0050     typedef typename odeint::unwrap_reference< Observer >::type observer_type;
0051
0052     stepper_type &st = stepper;
0053     observer_type &obs = observer;
0054     typedef typename unit_value_type<Time>::type time_type;
0055
0056     size_t steps = 0;
0057     Time current_dt = dt;
0058     while( true )
0059     {
0060         Time current_time = *start_time++;
0061         obs( start_state , current_time );
0062         if( start_time == end_time )
0063             break;
0064         while( less_with_sign( current_time , static_cast<time_type>(*start_time) , current_dt ) )
0065         {
0066             current_dt = min_abs( dt , *start_time - current_time );
0067             st.do_step( system , start_state , current_time , current_dt );
0068             current_time += current_dt;
0069             steps++;
0070         }
0071     }
0072     return steps;
0073 }
0074
0075 /*
0076  * integrate_times for controlled stepper
0077  */
0078 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
0079 size_t integrate_times(
0080         Stepper stepper , System system , State &start_state ,
0081         TimeIterator start_time , TimeIterator end_time , Time dt ,
0082         Observer observer , controlled_stepper_tag
0083 )
0084 {
0085     typename odeint::unwrap_reference< Observer >::type &obs = observer;
0086     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
0087     typedef typename unit_value_type<Time>::type time_type;
0088
0089     failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
0090     size_t steps = 0;
0091     while( true )
0092     {
0093         Time current_time = *start_time++;
0094         obs( start_state , current_time );
0095         if( start_time == end_time )
0096             break;
0097         while( less_with_sign( current_time , static_cast<time_type>(*start_time) , dt ) )
0098         {
0099             // adjust stepsize to end up exactly at the observation point
0100             Time current_dt = min_abs( dt , *start_time - current_time );
0101             if( st.try_step( system , start_state , current_time , current_dt ) == success )
0102             {
0103                 ++steps;
0104                 // successful step -> reset the fail counter, see #173
0105                 fail_checker.reset();
0106                 // continue with the original step size if dt was reduced due to observation
0107                 dt = max_abs( dt , current_dt );
0108             }
0109             else
0110             {
0111                 fail_checker();  // check for possible overflow of failed steps in step size adjustment
0112                 dt = current_dt;
0113             }
0114         }
0115     }
0116     return steps;
0117 }
0118
0119 /*
0120  * integrate_times for dense output stepper
0121  */
0122 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
0123 size_t integrate_times(
0124         Stepper stepper , System system , State &start_state ,
0125         TimeIterator start_time , TimeIterator end_time , Time dt ,
0126         Observer observer , dense_output_stepper_tag
0127 )
0128 {
0129     typename odeint::unwrap_reference< Observer >::type &obs = observer;
0130     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
0131
0132     typedef typename unit_value_type<Time>::type time_type;
0133
0134     if( start_time == end_time )
0135         return 0;
0136
0137     TimeIterator last_time_iterator = end_time;
0138     --last_time_iterator;
0139     Time last_time_point = static_cast<time_type>(*last_time_iterator);
0140
0141     st.initialize( start_state , *start_time , dt );
0142     obs( start_state , *start_time++ );
0143
0144     size_t count = 0;
0145     while( start_time != end_time )
0146     {
0147         while( ( start_time != end_time ) && less_eq_with_sign( static_cast<time_type>(*start_time) , st.current_time() , st.current_time_step() ) )
0148         {
0149             st.calc_state( *start_time , start_state );
0150             obs( start_state , *start_time );
0151             start_time++;
0152         }
0153
0154         // we have not reached the end, do another real step
0155         if( less_eq_with_sign( st.current_time() + st.current_time_step() ,
0156                                last_time_point ,
0157                                st.current_time_step() ) )
0158         {
0159             st.do_step( system );
0160             ++count;
0161         }
0162         else if( start_time != end_time )
0163         { // do the last step ending exactly on the end point
0164             st.initialize( st.current_state() , st.current_time() , last_time_point - st.current_time() );
0165             st.do_step( system );
0166             ++count;
0167         }
0168     }
0169     return count;
0170 }
0171
0172
0173 } // namespace detail
0174 } // namespace odeint
0175 } // namespace numeric
0176 } // namespace boost
0177
0178
0179 #endif // BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED