3 boost/numeric/odeint/integrate/detail/integrate_times.hpp
6 Default integrate times implementation.
9 Copyright 2011-2015 Mario Mulansky
10 Copyright 2012 Karsten Ahnert
11 Copyright 2012 Christoph Koke
13 Distributed under the Boost Software License, Version 1.0.
14 (See accompanying file LICENSE_1_0.txt or
15 copy at http://www.boost.org/LICENSE_1_0.txt)
19 #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
20 #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
24 #include <boost/config.hpp>
25 #include <boost/throw_exception.hpp>
26 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
27 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
28 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
29 #include <boost/numeric/odeint/integrate/max_step_checker.hpp>
40 * integrate_times for simple stepper
42 template<class Stepper, class System, class State, class TimeIterator, class Time, class Observer>
43 size_t integrate_times(
44 Stepper stepper , System system , State &start_state ,
45 TimeIterator start_time , TimeIterator end_time , Time dt ,
46 Observer observer , stepper_tag
49 typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
50 typedef typename odeint::unwrap_reference< Observer >::type observer_type;
52 stepper_type &st = stepper;
53 observer_type &obs = observer;
54 typedef typename unit_value_type<Time>::type time_type;
60 Time current_time = *start_time++;
61 obs( start_state , current_time );
62 if( start_time == end_time )
64 while( less_with_sign( current_time , static_cast<time_type>(*start_time) , current_dt ) )
66 current_dt = min_abs( dt , *start_time - current_time );
67 st.do_step( system , start_state , current_time , current_dt );
68 current_time += current_dt;
76 * integrate_times for controlled stepper
78 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
79 size_t integrate_times(
80 Stepper stepper , System system , State &start_state ,
81 TimeIterator start_time , TimeIterator end_time , Time dt ,
82 Observer observer , controlled_stepper_tag
85 typename odeint::unwrap_reference< Observer >::type &obs = observer;
86 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
87 typedef typename unit_value_type<Time>::type time_type;
89 failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails
93 Time current_time = *start_time++;
94 obs( start_state , current_time );
95 if( start_time == end_time )
97 while( less_with_sign( current_time , static_cast<time_type>(*start_time) , dt ) )
99 // adjust stepsize to end up exactly at the observation point
100 Time current_dt = min_abs( dt , *start_time - current_time );
101 if( st.try_step( system , start_state , current_time , current_dt ) == success )
104 // successful step -> reset the fail counter, see #173
105 fail_checker.reset();
106 // continue with the original step size if dt was reduced due to observation
107 dt = max_abs( dt , current_dt );
111 fail_checker(); // check for possible overflow of failed steps in step size adjustment
120 * integrate_times for dense output stepper
122 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
123 size_t integrate_times(
124 Stepper stepper , System system , State &start_state ,
125 TimeIterator start_time , TimeIterator end_time , Time dt ,
126 Observer observer , dense_output_stepper_tag
129 typename odeint::unwrap_reference< Observer >::type &obs = observer;
130 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
132 typedef typename unit_value_type<Time>::type time_type;
134 if( start_time == end_time )
137 TimeIterator last_time_iterator = end_time;
138 --last_time_iterator;
139 Time last_time_point = static_cast<time_type>(*last_time_iterator);
141 st.initialize( start_state , *start_time , dt );
142 obs( start_state , *start_time++ );
145 while( start_time != end_time )
147 while( ( start_time != end_time ) && less_eq_with_sign( static_cast<time_type>(*start_time) , st.current_time() , st.current_time_step() ) )
149 st.calc_state( *start_time , start_state );
150 obs( start_state , *start_time );
154 // we have not reached the end, do another real step
155 if( less_eq_with_sign( st.current_time() + st.current_time_step() ,
157 st.current_time_step() ) )
159 st.do_step( system );
162 else if( start_time != end_time )
163 { // do the last step ending exactly on the end point
164 st.initialize( st.current_state() , st.current_time() , last_time_point - st.current_time() );
165 st.do_step( system );
173 } // namespace detail
174 } // namespace odeint
175 } // namespace numeric
179 #endif // BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED