]>
Commit | Line | Data |
---|---|---|
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> | |
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 > | |
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 > | |
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 |