Chris@16
|
1 /*
|
Chris@16
|
2 [auto_generated]
|
Chris@16
|
3 boost/numeric/odeint/integrate/detail/integrate_times.hpp
|
Chris@16
|
4
|
Chris@16
|
5 [begin_description]
|
Chris@16
|
6 Default integrate times implementation.
|
Chris@16
|
7 [end_description]
|
Chris@16
|
8
|
Chris@101
|
9 Copyright 2011-2012 Mario Mulansky
|
Chris@101
|
10 Copyright 2012 Karsten Ahnert
|
Chris@101
|
11 Copyright 2012 Christoph Koke
|
Chris@16
|
12
|
Chris@16
|
13 Distributed under the Boost Software License, Version 1.0.
|
Chris@16
|
14 (See accompanying file LICENSE_1_0.txt or
|
Chris@16
|
15 copy at http://www.boost.org/LICENSE_1_0.txt)
|
Chris@16
|
16 */
|
Chris@16
|
17
|
Chris@16
|
18
|
Chris@16
|
19 #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
|
Chris@16
|
20 #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
|
Chris@16
|
21
|
Chris@16
|
22 #include <stdexcept>
|
Chris@16
|
23
|
Chris@16
|
24 #include <boost/config.hpp>
|
Chris@101
|
25 #include <boost/throw_exception.hpp>
|
Chris@16
|
26 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
|
Chris@16
|
27 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
|
Chris@16
|
28 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
|
Chris@16
|
29
|
Chris@16
|
30
|
Chris@16
|
31 namespace boost {
|
Chris@16
|
32 namespace numeric {
|
Chris@16
|
33 namespace odeint {
|
Chris@16
|
34 namespace detail {
|
Chris@16
|
35
|
Chris@16
|
36
|
Chris@16
|
37
|
Chris@16
|
38 /*
|
Chris@16
|
39 * integrate_times for simple stepper
|
Chris@16
|
40 */
|
Chris@16
|
41 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
|
Chris@16
|
42 size_t integrate_times(
|
Chris@16
|
43 Stepper stepper , System system , State &start_state ,
|
Chris@16
|
44 TimeIterator start_time , TimeIterator end_time , Time dt ,
|
Chris@16
|
45 Observer observer , stepper_tag
|
Chris@16
|
46 )
|
Chris@16
|
47 {
|
Chris@16
|
48 typename odeint::unwrap_reference< Observer >::type &obs = observer;
|
Chris@101
|
49 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
|
Chris@101
|
50 typedef typename unit_value_type<Time>::type time_type;
|
Chris@16
|
51
|
Chris@16
|
52 size_t steps = 0;
|
Chris@16
|
53 Time current_dt = dt;
|
Chris@16
|
54 while( true )
|
Chris@16
|
55 {
|
Chris@16
|
56 Time current_time = *start_time++;
|
Chris@16
|
57 obs( start_state , current_time );
|
Chris@16
|
58 if( start_time == end_time )
|
Chris@16
|
59 break;
|
Chris@101
|
60 while( less_with_sign( current_time , static_cast<time_type>(*start_time) , current_dt ) )
|
Chris@16
|
61 {
|
Chris@101
|
62 current_dt = min_abs( dt , *start_time - current_time );
|
Chris@101
|
63 st.do_step( system , start_state , current_time , current_dt );
|
Chris@16
|
64 current_time += current_dt;
|
Chris@16
|
65 steps++;
|
Chris@16
|
66 }
|
Chris@16
|
67 }
|
Chris@16
|
68 return steps;
|
Chris@16
|
69 }
|
Chris@16
|
70
|
Chris@16
|
71 /*
|
Chris@16
|
72 * integrate_times for controlled stepper
|
Chris@16
|
73 */
|
Chris@16
|
74 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
|
Chris@16
|
75 size_t integrate_times(
|
Chris@16
|
76 Stepper stepper , System system , State &start_state ,
|
Chris@16
|
77 TimeIterator start_time , TimeIterator end_time , Time dt ,
|
Chris@16
|
78 Observer observer , controlled_stepper_tag
|
Chris@16
|
79 )
|
Chris@16
|
80 {
|
Chris@16
|
81 typename odeint::unwrap_reference< Observer >::type &obs = observer;
|
Chris@101
|
82 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
|
Chris@101
|
83 typedef typename unit_value_type<Time>::type time_type;
|
Chris@16
|
84
|
Chris@16
|
85 const size_t max_attempts = 1000;
|
Chris@16
|
86 const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found.";
|
Chris@16
|
87 size_t steps = 0;
|
Chris@16
|
88 while( true )
|
Chris@16
|
89 {
|
Chris@16
|
90 size_t fail_steps = 0;
|
Chris@16
|
91 Time current_time = *start_time++;
|
Chris@16
|
92 obs( start_state , current_time );
|
Chris@16
|
93 if( start_time == end_time )
|
Chris@16
|
94 break;
|
Chris@101
|
95 while( less_with_sign( current_time , static_cast<time_type>(*start_time) , dt ) )
|
Chris@16
|
96 {
|
Chris@101
|
97 // adjust stepsize to end up exactly at the observation point
|
Chris@101
|
98 Time current_dt = min_abs( dt , *start_time - current_time );
|
Chris@101
|
99 if( st.try_step( system , start_state , current_time , current_dt ) == success )
|
Chris@16
|
100 {
|
Chris@16
|
101 ++steps;
|
Chris@101
|
102 // continue with the original step size if dt was reduced due to observation
|
Chris@101
|
103 dt = max_abs( dt , current_dt );
|
Chris@16
|
104 }
|
Chris@16
|
105 else
|
Chris@16
|
106 {
|
Chris@16
|
107 ++fail_steps;
|
Chris@101
|
108 dt = current_dt;
|
Chris@16
|
109 }
|
Chris@101
|
110 if( fail_steps == max_attempts ) BOOST_THROW_EXCEPTION( std::overflow_error( error_string ));
|
Chris@16
|
111 }
|
Chris@16
|
112 }
|
Chris@16
|
113 return steps;
|
Chris@16
|
114 }
|
Chris@16
|
115
|
Chris@16
|
116 /*
|
Chris@16
|
117 * integrate_times for dense output stepper
|
Chris@16
|
118 */
|
Chris@16
|
119 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
|
Chris@16
|
120 size_t integrate_times(
|
Chris@16
|
121 Stepper stepper , System system , State &start_state ,
|
Chris@16
|
122 TimeIterator start_time , TimeIterator end_time , Time dt ,
|
Chris@16
|
123 Observer observer , dense_output_stepper_tag
|
Chris@16
|
124 )
|
Chris@16
|
125 {
|
Chris@16
|
126 typename odeint::unwrap_reference< Observer >::type &obs = observer;
|
Chris@101
|
127 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
|
Chris@101
|
128 typedef typename unit_value_type<Time>::type time_type;
|
Chris@16
|
129
|
Chris@16
|
130 if( start_time == end_time )
|
Chris@16
|
131 return 0;
|
Chris@16
|
132
|
Chris@101
|
133 Time last_time_point = static_cast<time_type>(*(end_time-1));
|
Chris@16
|
134
|
Chris@101
|
135 st.initialize( start_state , *start_time , dt );
|
Chris@16
|
136 obs( start_state , *start_time++ );
|
Chris@16
|
137
|
Chris@16
|
138 size_t count = 0;
|
Chris@16
|
139 while( start_time != end_time )
|
Chris@16
|
140 {
|
Chris@101
|
141 while( ( start_time != end_time ) && less_eq_with_sign( static_cast<time_type>(*start_time) , st.current_time() , st.current_time_step() ) )
|
Chris@16
|
142 {
|
Chris@101
|
143 st.calc_state( *start_time , start_state );
|
Chris@16
|
144 obs( start_state , *start_time );
|
Chris@16
|
145 start_time++;
|
Chris@16
|
146 }
|
Chris@16
|
147
|
Chris@16
|
148 // we have not reached the end, do another real step
|
Chris@101
|
149 if( less_eq_with_sign( st.current_time() + st.current_time_step() ,
|
Chris@16
|
150 last_time_point ,
|
Chris@101
|
151 st.current_time_step() ) )
|
Chris@16
|
152 {
|
Chris@101
|
153 st.do_step( system );
|
Chris@16
|
154 ++count;
|
Chris@16
|
155 }
|
Chris@16
|
156 else if( start_time != end_time )
|
Chris@16
|
157 { // do the last step ending exactly on the end point
|
Chris@101
|
158 st.initialize( st.current_state() , st.current_time() , last_time_point - st.current_time() );
|
Chris@101
|
159 st.do_step( system );
|
Chris@16
|
160 ++count;
|
Chris@16
|
161 }
|
Chris@16
|
162 }
|
Chris@16
|
163 return count;
|
Chris@16
|
164 }
|
Chris@16
|
165
|
Chris@16
|
166
|
Chris@16
|
167 } // namespace detail
|
Chris@16
|
168 } // namespace odeint
|
Chris@16
|
169 } // namespace numeric
|
Chris@16
|
170 } // namespace boost
|
Chris@16
|
171
|
Chris@16
|
172
|
Chris@16
|
173 #endif // BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED
|