1 /*
2 [auto_generated]
3 boost/numeric/odeint/integrate/detail/integrate_times.hpp
4
5 [begin_description]
6 Default integrate times implementation.
7 [end_description]
8
9 Copyright 2011-2015 Mario Mulansky
10 Copyright 2012 Karsten Ahnert
11 Copyright 2012 Christoph Koke
12
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)
16 */
17
18
19 #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
20 #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
21
22 #include <stdexcept>
23
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>
30
31
32 namespace boost {
33 namespace numeric {
34 namespace odeint {
35 namespace detail {
36
37
38
39 /*
40 * integrate_times for simple stepper
41 */
42 template<class Stepper, class System, class State, class TimeIterator, class Time, class Observer>
integrate_times(Stepper stepper,System system,State & start_state,TimeIterator start_time,TimeIterator end_time,Time dt,Observer observer,stepper_tag)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
47 )
48 {
49 typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
50 typedef typename odeint::unwrap_reference< Observer >::type observer_type;
51
52 stepper_type &st = stepper;
53 observer_type &obs = observer;
54 typedef typename unit_value_type<Time>::type time_type;
55
56 size_t steps = 0;
57 Time current_dt = dt;
58 while( true )
59 {
60 Time current_time = *start_time++;
61 obs( start_state , current_time );
62 if( start_time == end_time )
63 break;
64 while( less_with_sign( current_time , static_cast<time_type>(*start_time) , current_dt ) )
65 {
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;
69 steps++;
70 }
71 }
72 return steps;
73 }
74
75 /*
76 * integrate_times for controlled stepper
77 */
78 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
integrate_times(Stepper stepper,System system,State & start_state,TimeIterator start_time,TimeIterator end_time,Time dt,Observer observer,controlled_stepper_tag)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
83 )
84 {
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;
88
89 failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails
90 size_t steps = 0;
91 while( true )
92 {
93 Time current_time = *start_time++;
94 obs( start_state , current_time );
95 if( start_time == end_time )
96 break;
97 while( less_with_sign( current_time , static_cast<time_type>(*start_time) , dt ) )
98 {
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 )
102 {
103 ++steps;
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 );
108 }
109 else
110 {
111 fail_checker(); // check for possible overflow of failed steps in step size adjustment
112 dt = current_dt;
113 }
114 }
115 }
116 return steps;
117 }
118
119 /*
120 * integrate_times for dense output stepper
121 */
122 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
integrate_times(Stepper stepper,System system,State & start_state,TimeIterator start_time,TimeIterator end_time,Time dt,Observer observer,dense_output_stepper_tag)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
127 )
128 {
129 typename odeint::unwrap_reference< Observer >::type &obs = observer;
130 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
131
132 typedef typename unit_value_type<Time>::type time_type;
133
134 if( start_time == end_time )
135 return 0;
136
137 TimeIterator last_time_iterator = end_time;
138 --last_time_iterator;
139 Time last_time_point = static_cast<time_type>(*last_time_iterator);
140
141 st.initialize( start_state , *start_time , dt );
142 obs( start_state , *start_time++ );
143
144 size_t count = 0;
145 while( start_time != end_time )
146 {
147 while( ( start_time != end_time ) && less_eq_with_sign( static_cast<time_type>(*start_time) , st.current_time() , st.current_time_step() ) )
148 {
149 st.calc_state( *start_time , start_state );
150 obs( start_state , *start_time );
151 start_time++;
152 }
153
154 // we have not reached the end, do another real step
155 if( less_eq_with_sign( st.current_time() + st.current_time_step() ,
156 last_time_point ,
157 st.current_time_step() ) )
158 {
159 st.do_step( system );
160 ++count;
161 }
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 );
166 ++count;
167 }
168 }
169 return count;
170 }
171
172
173 } // namespace detail
174 } // namespace odeint
175 } // namespace numeric
176 } // namespace boost
177
178
179 #endif // BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED
180