]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | /* |
2 | [auto_generated] | |
3 | boost/numeric/odeint/integrate/detail/integrate_const.hpp | |
4 | ||
5 | [begin_description] | |
6 | integrate const implementation | |
7 | [end_description] | |
8 | ||
9 | Copyright 2012-2015 Mario Mulansky | |
10 | Copyright 2012 Christoph Koke | |
11 | Copyright 2012 Karsten Ahnert | |
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 | #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED | |
19 | #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED | |
20 | ||
21 | #include <boost/numeric/odeint/util/unwrap_reference.hpp> | |
22 | #include <boost/numeric/odeint/stepper/stepper_categories.hpp> | |
23 | #include <boost/numeric/odeint/util/unit_helper.hpp> | |
24 | #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp> | |
25 | ||
26 | #include <boost/numeric/odeint/util/detail/less_with_sign.hpp> | |
27 | ||
28 | namespace boost { | |
29 | namespace numeric { | |
30 | namespace odeint { | |
31 | namespace detail { | |
32 | ||
33 | // forward declaration | |
34 | template< class Stepper , class System , class State , class Time , class Observer > | |
35 | size_t integrate_adaptive( | |
36 | Stepper stepper , System system , State &start_state , | |
37 | Time &start_time , Time end_time , Time &dt , | |
38 | Observer observer , controlled_stepper_tag | |
39 | ); | |
40 | ||
41 | ||
42 | template< class Stepper , class System , class State , class Time , class Observer > | |
43 | size_t integrate_const( | |
44 | Stepper stepper , System system , State &start_state , | |
45 | Time start_time , Time end_time , Time dt , | |
46 | Observer observer , stepper_tag | |
47 | ) | |
48 | { | |
49 | ||
50 | typename odeint::unwrap_reference< Observer >::type &obs = observer; | |
51 | typename odeint::unwrap_reference< Stepper >::type &st = stepper; | |
52 | ||
53 | Time time = start_time; | |
54 | int step = 0; | |
55 | // cast time+dt explicitely in case of expression templates (e.g. multiprecision) | |
56 | while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) ) | |
57 | { | |
58 | obs( start_state , time ); | |
59 | st.do_step( system , start_state , time , dt ); | |
60 | // direct computation of the time avoids error propagation happening when using time += dt | |
61 | // we need clumsy type analysis to get boost units working here | |
62 | ++step; | |
63 | time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * dt; | |
64 | } | |
65 | obs( start_state , time ); | |
66 | ||
67 | return step; | |
68 | } | |
69 | ||
70 | ||
71 | ||
72 | template< class Stepper , class System , class State , class Time , class Observer > | |
73 | size_t integrate_const( | |
74 | Stepper stepper , System system , State &start_state , | |
75 | Time start_time , Time end_time , Time dt , | |
76 | Observer observer , controlled_stepper_tag | |
77 | ) | |
78 | { | |
79 | typename odeint::unwrap_reference< Observer >::type &obs = observer; | |
80 | ||
81 | Time time = start_time; | |
82 | const Time time_step = dt; | |
83 | int real_steps = 0; | |
84 | int step = 0; | |
85 | ||
86 | while( less_eq_with_sign( static_cast<Time>(time+time_step) , end_time , dt ) ) | |
87 | { | |
88 | obs( start_state , time ); | |
89 | // integrate_adaptive_checked uses the given checker to throw if an overflow occurs | |
90 | real_steps += detail::integrate_adaptive(stepper, system, start_state, time, | |
91 | static_cast<Time>(time + time_step), dt, | |
92 | null_observer(), controlled_stepper_tag()); | |
93 | // direct computation of the time avoids error propagation happening when using time += dt | |
94 | // we need clumsy type analysis to get boost units working here | |
95 | step++; | |
96 | time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * time_step; | |
97 | } | |
98 | obs( start_state , time ); | |
99 | ||
100 | return real_steps; | |
101 | } | |
102 | ||
103 | ||
104 | template< class Stepper , class System , class State , class Time , class Observer > | |
105 | size_t integrate_const( | |
106 | Stepper stepper , System system , State &start_state , | |
107 | Time start_time , Time end_time , Time dt , | |
108 | Observer observer , dense_output_stepper_tag | |
109 | ) | |
110 | { | |
111 | typename odeint::unwrap_reference< Observer >::type &obs = observer; | |
112 | typename odeint::unwrap_reference< Stepper >::type &st = stepper; | |
113 | ||
114 | Time time = start_time; | |
115 | ||
116 | st.initialize( start_state , time , dt ); | |
117 | obs( start_state , time ); | |
118 | time += dt; | |
119 | ||
120 | int obs_step( 1 ); | |
121 | int real_step( 0 ); | |
122 | ||
123 | while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) ) | |
124 | { | |
125 | while( less_eq_with_sign( time , st.current_time() , dt ) ) | |
126 | { | |
127 | st.calc_state( time , start_state ); | |
128 | obs( start_state , time ); | |
129 | ++obs_step; | |
130 | // direct computation of the time avoids error propagation happening when using time += dt | |
131 | // we need clumsy type analysis to get boost units working here | |
132 | time = start_time + static_cast< typename unit_value_type<Time>::type >(obs_step) * dt; | |
133 | } | |
134 | // we have not reached the end, do another real step | |
135 | if( less_with_sign( static_cast<Time>(st.current_time()+st.current_time_step()) , | |
136 | end_time , | |
137 | st.current_time_step() ) ) | |
138 | { | |
139 | while( less_eq_with_sign( st.current_time() , time , dt ) ) | |
140 | { | |
141 | st.do_step( system ); | |
142 | ++real_step; | |
143 | } | |
144 | } | |
145 | else if( less_with_sign( st.current_time() , end_time , st.current_time_step() ) ) | |
146 | { // do the last step ending exactly on the end point | |
147 | st.initialize( st.current_state() , st.current_time() , end_time - st.current_time() ); | |
148 | st.do_step( system ); | |
149 | ++real_step; | |
150 | } | |
151 | ||
152 | } | |
153 | // last observation, if we are still in observation interval | |
154 | // might happen due to finite precision problems when computing the the time | |
155 | if( less_eq_with_sign( time , end_time , dt ) ) | |
156 | { | |
157 | st.calc_state( time , start_state ); | |
158 | obs( start_state , time ); | |
159 | } | |
160 | ||
161 | return real_step; | |
162 | } | |
163 | ||
164 | ||
165 | } } } } | |
166 | ||
167 | #endif |