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-2012 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
30
31 namespace boost {
32 namespace numeric {
33 namespace odeint {
34 namespace detail {
35
36
37
38 /*
39 * integrate_times for simple stepper
40 */
41 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)42 size_t integrate_times(
43 Stepper stepper , System system , State &start_state ,
44 TimeIterator start_time , TimeIterator end_time , Time dt ,
45 Observer observer , stepper_tag
46 )
47 {
48 typename odeint::unwrap_reference< Observer >::type &obs = observer;
49 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
50 typedef typename unit_value_type<Time>::type time_type;
51
52 size_t steps = 0;
53 Time current_dt = dt;
54 while( true )
55 {
56 Time current_time = *start_time++;
57 obs( start_state , current_time );
58 if( start_time == end_time )
59 break;
60 while( less_with_sign( current_time , static_cast<time_type>(*start_time) , current_dt ) )
61 {
62 current_dt = min_abs( dt , *start_time - current_time );
63 st.do_step( system , start_state , current_time , current_dt );
64 current_time += current_dt;
65 steps++;
66 }
67 }
68 return steps;
69 }
70
71 /*
72 * integrate_times for controlled stepper
73 */
74 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)75 size_t integrate_times(
76 Stepper stepper , System system , State &start_state ,
77 TimeIterator start_time , TimeIterator end_time , Time dt ,
78 Observer observer , controlled_stepper_tag
79 )
80 {
81 typename odeint::unwrap_reference< Observer >::type &obs = observer;
82 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
83 typedef typename unit_value_type<Time>::type time_type;
84
85 const size_t max_attempts = 1000;
86 const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found.";
87 size_t steps = 0;
88 while( true )
89 {
90 size_t fail_steps = 0;
91 Time current_time = *start_time++;
92 obs( start_state , current_time );
93 if( start_time == end_time )
94 break;
95 while( less_with_sign( current_time , static_cast<time_type>(*start_time) , dt ) )
96 {
97 // adjust stepsize to end up exactly at the observation point
98 Time current_dt = min_abs( dt , *start_time - current_time );
99 if( st.try_step( system , start_state , current_time , current_dt ) == success )
100 {
101 ++steps;
102 // continue with the original step size if dt was reduced due to observation
103 dt = max_abs( dt , current_dt );
104 }
105 else
106 {
107 ++fail_steps;
108 dt = current_dt;
109 }
110 if( fail_steps == max_attempts ) BOOST_THROW_EXCEPTION( std::overflow_error( error_string ));
111 }
112 }
113 return steps;
114 }
115
116 /*
117 * integrate_times for dense output stepper
118 */
119 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)120 size_t integrate_times(
121 Stepper stepper , System system , State &start_state ,
122 TimeIterator start_time , TimeIterator end_time , Time dt ,
123 Observer observer , dense_output_stepper_tag
124 )
125 {
126 typename odeint::unwrap_reference< Observer >::type &obs = observer;
127 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
128 typedef typename unit_value_type<Time>::type time_type;
129
130 if( start_time == end_time )
131 return 0;
132
133 TimeIterator last_time_iterator = end_time;
134 --last_time_iterator;
135 Time last_time_point = static_cast<time_type>(*last_time_iterator);
136
137 st.initialize( start_state , *start_time , dt );
138 obs( start_state , *start_time++ );
139
140 size_t count = 0;
141 while( start_time != end_time )
142 {
143 while( ( start_time != end_time ) && less_eq_with_sign( static_cast<time_type>(*start_time) , st.current_time() , st.current_time_step() ) )
144 {
145 st.calc_state( *start_time , start_state );
146 obs( start_state , *start_time );
147 start_time++;
148 }
149
150 // we have not reached the end, do another real step
151 if( less_eq_with_sign( st.current_time() + st.current_time_step() ,
152 last_time_point ,
153 st.current_time_step() ) )
154 {
155 st.do_step( system );
156 ++count;
157 }
158 else if( start_time != end_time )
159 { // do the last step ending exactly on the end point
160 st.initialize( st.current_state() , st.current_time() , last_time_point - st.current_time() );
161 st.do_step( system );
162 ++count;
163 }
164 }
165 return count;
166 }
167
168
169 } // namespace detail
170 } // namespace odeint
171 } // namespace numeric
172 } // namespace boost
173
174
175 #endif // BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED
176