]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | /* |
2 | [auto_generated] | |
3 | boost/numeric/odeint/stepper/runge_kutta_fehlberg87.hpp | |
4 | ||
5 | [begin_description] | |
6 | Implementation of the Runge-Kutta-Fehlberg stepper with the generic stepper. | |
7 | [end_description] | |
8 | ||
9 | Copyright 2011-2013 Mario Mulansky | |
10 | Copyright 2012-2013 Karsten Ahnert | |
11 | ||
12 | Distributed under the Boost Software License, Version 1.0. | |
13 | (See accompanying file LICENSE_1_0.txt or | |
14 | copy at http://www.boost.org/LICENSE_1_0.txt) | |
15 | */ | |
16 | ||
17 | ||
18 | #ifndef BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_FEHLBERG87_HPP_INCLUDED | |
19 | #define BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_FEHLBERG87_HPP_INCLUDED | |
20 | ||
21 | ||
22 | #include <boost/fusion/container/vector.hpp> | |
23 | #include <boost/fusion/container/generation/make_vector.hpp> | |
24 | ||
25 | #include <boost/numeric/odeint/stepper/explicit_error_generic_rk.hpp> | |
26 | #include <boost/numeric/odeint/algebra/range_algebra.hpp> | |
27 | #include <boost/numeric/odeint/algebra/default_operations.hpp> | |
28 | #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> | |
29 | #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> | |
30 | ||
31 | #include <boost/array.hpp> | |
32 | ||
33 | #include <boost/numeric/odeint/util/state_wrapper.hpp> | |
34 | #include <boost/numeric/odeint/util/is_resizeable.hpp> | |
35 | #include <boost/numeric/odeint/util/resizer.hpp> | |
36 | ||
37 | ||
38 | ||
39 | ||
40 | namespace boost { | |
41 | namespace numeric { | |
42 | namespace odeint { | |
43 | ||
44 | ||
45 | #ifndef DOXYGEN_SKIP | |
46 | template< class Value = double > | |
47 | struct rk78_coefficients_a1 : boost::array< Value , 1 > | |
48 | { | |
49 | rk78_coefficients_a1( void ) | |
50 | { | |
51 | (*this)[0] = static_cast< Value >( 2 )/static_cast< Value >( 27 ); | |
52 | } | |
53 | }; | |
54 | ||
55 | template< class Value = double > | |
56 | struct rk78_coefficients_a2 : boost::array< Value , 2 > | |
57 | { | |
58 | rk78_coefficients_a2( void ) | |
59 | { | |
60 | (*this)[0] = static_cast< Value >( 1 )/static_cast< Value >( 36 ); | |
61 | (*this)[1] = static_cast< Value >( 1 )/static_cast< Value >( 12 ); | |
62 | } | |
63 | }; | |
64 | ||
65 | ||
66 | template< class Value = double > | |
67 | struct rk78_coefficients_a3 : boost::array< Value , 3 > | |
68 | { | |
69 | rk78_coefficients_a3( void ) | |
70 | { | |
71 | (*this)[0] = static_cast< Value >( 1 )/static_cast< Value >( 24 ); | |
72 | (*this)[1] = static_cast< Value >( 0 ); | |
73 | (*this)[2] = static_cast< Value >( 1 )/static_cast< Value >( 8 ); | |
74 | } | |
75 | }; | |
76 | ||
77 | template< class Value = double > | |
78 | struct rk78_coefficients_a4 : boost::array< Value , 4 > | |
79 | { | |
80 | rk78_coefficients_a4( void ) | |
81 | { | |
82 | (*this)[0] = static_cast< Value >( 5 )/static_cast< Value >( 12 ); | |
83 | (*this)[1] = static_cast< Value >( 0 ); | |
84 | (*this)[2] = static_cast< Value >( -25 )/static_cast< Value >( 16 ); | |
85 | (*this)[3] = static_cast< Value >( 25 )/static_cast< Value >( 16 ); | |
86 | } | |
87 | }; | |
88 | ||
89 | template< class Value = double > | |
90 | struct rk78_coefficients_a5 : boost::array< Value , 5 > | |
91 | { | |
92 | rk78_coefficients_a5( void ) | |
93 | { | |
94 | (*this)[0] = static_cast< Value >( 1 )/static_cast< Value >( 20 ); | |
95 | (*this)[1] = static_cast< Value >( 0 ); | |
96 | (*this)[2] = static_cast< Value >( 0 ); | |
97 | (*this)[3] = static_cast< Value >( 1 )/static_cast< Value >( 4 ); | |
98 | (*this)[4] = static_cast< Value >( 1 )/static_cast< Value >( 5 ); | |
99 | } | |
100 | }; | |
101 | ||
102 | ||
103 | template< class Value = double > | |
104 | struct rk78_coefficients_a6 : boost::array< Value , 6 > | |
105 | { | |
106 | rk78_coefficients_a6( void ) | |
107 | { | |
108 | (*this)[0] = static_cast< Value >( -25 )/static_cast< Value >( 108 ); | |
109 | (*this)[1] = static_cast< Value >( 0 ); | |
110 | (*this)[2] = static_cast< Value >( 0 ); | |
111 | (*this)[3] = static_cast< Value >( 125 )/static_cast< Value >( 108 ); | |
112 | (*this)[4] = static_cast< Value >( -65 )/static_cast< Value >( 27 ); | |
113 | (*this)[5] = static_cast< Value >( 125 )/static_cast< Value >( 54 ); | |
114 | } | |
115 | }; | |
116 | ||
117 | template< class Value = double > | |
118 | struct rk78_coefficients_a7 : boost::array< Value , 7 > | |
119 | { | |
120 | rk78_coefficients_a7( void ) | |
121 | { | |
122 | (*this)[0] = static_cast< Value >( 31 )/static_cast< Value >( 300 ); | |
123 | (*this)[1] = static_cast< Value >( 0 ); | |
124 | (*this)[2] = static_cast< Value >( 0 ); | |
125 | (*this)[3] = static_cast< Value >( 0 ); | |
126 | (*this)[4] = static_cast< Value >( 61 )/static_cast< Value >( 225 ); | |
127 | (*this)[5] = static_cast< Value >( -2 )/static_cast< Value >( 9 ); | |
128 | (*this)[6] = static_cast< Value >( 13 )/static_cast< Value >( 900 ); | |
129 | } | |
130 | }; | |
131 | ||
132 | template< class Value = double > | |
133 | struct rk78_coefficients_a8 : boost::array< Value , 8 > | |
134 | { | |
135 | rk78_coefficients_a8( void ) | |
136 | { | |
137 | (*this)[0] = static_cast< Value >( 2 ); | |
138 | (*this)[1] = static_cast< Value >( 0 ); | |
139 | (*this)[2] = static_cast< Value >( 0 ); | |
140 | (*this)[3] = static_cast< Value >( -53 )/static_cast< Value >( 6 ); | |
141 | (*this)[4] = static_cast< Value >( 704 )/static_cast< Value >( 45 ); | |
142 | (*this)[5] = static_cast< Value >( -107 )/static_cast< Value >( 9 ); | |
143 | (*this)[6] = static_cast< Value >( 67 )/static_cast< Value >( 90 ); | |
144 | (*this)[7] = static_cast< Value >( 3 ); | |
145 | } | |
146 | }; | |
147 | ||
148 | template< class Value = double > | |
149 | struct rk78_coefficients_a9 : boost::array< Value , 9 > | |
150 | { | |
151 | rk78_coefficients_a9( void ) | |
152 | { | |
153 | (*this)[0] = static_cast< Value >( -91 )/static_cast< Value >( 108 ); | |
154 | (*this)[1] = static_cast< Value >( 0 ); | |
155 | (*this)[2] = static_cast< Value >( 0 ); | |
156 | (*this)[3] = static_cast< Value >( 23 )/static_cast< Value >( 108 ); | |
157 | (*this)[4] = static_cast< Value >( -976 )/static_cast< Value >( 135 ); | |
158 | (*this)[5] = static_cast< Value >( 311 )/static_cast< Value >( 54 ); | |
159 | (*this)[6] = static_cast< Value >( -19 )/static_cast< Value >( 60 ); | |
160 | (*this)[7] = static_cast< Value >( 17 )/static_cast< Value >( 6 ); | |
161 | (*this)[8] = static_cast< Value >( -1 )/static_cast< Value >( 12 ); | |
162 | } | |
163 | }; | |
164 | ||
165 | template< class Value = double > | |
166 | struct rk78_coefficients_a10 : boost::array< Value , 10 > | |
167 | { | |
168 | rk78_coefficients_a10( void ) | |
169 | { | |
170 | (*this)[0] = static_cast< Value >( 2383 )/static_cast< Value >( 4100 ); | |
171 | (*this)[1] = static_cast< Value >( 0 ); | |
172 | (*this)[2] = static_cast< Value >( 0 ); | |
173 | (*this)[3] = static_cast< Value >( -341 )/static_cast< Value >( 164 ); | |
174 | (*this)[4] = static_cast< Value >( 4496 )/static_cast< Value >( 1025 ); | |
175 | (*this)[5] = static_cast< Value >( -301 )/static_cast< Value >( 82 ); | |
176 | (*this)[6] = static_cast< Value >( 2133 )/static_cast< Value >( 4100 ); | |
177 | (*this)[7] = static_cast< Value >( 45 )/static_cast< Value >( 82 ); | |
178 | (*this)[8] = static_cast< Value >( 45 )/static_cast< Value >( 164 ); | |
179 | (*this)[9] = static_cast< Value >( 18 )/static_cast< Value >( 41 ); | |
180 | } | |
181 | }; | |
182 | ||
183 | template< class Value = double > | |
184 | struct rk78_coefficients_a11 : boost::array< Value , 11 > | |
185 | { | |
186 | rk78_coefficients_a11( void ) | |
187 | { | |
188 | (*this)[0] = static_cast< Value >( 3 )/static_cast< Value >( 205 ); | |
189 | (*this)[1] = static_cast< Value >( 0 ); | |
190 | (*this)[2] = static_cast< Value >( 0 ); | |
191 | (*this)[3] = static_cast< Value >( 0 ); | |
192 | (*this)[4] = static_cast< Value >( 0 ); | |
193 | (*this)[5] = static_cast< Value >( -6 )/static_cast< Value >( 41 ); | |
194 | (*this)[6] = static_cast< Value >( -3 )/static_cast< Value >( 205 ); | |
195 | (*this)[7] = static_cast< Value >( -3 )/static_cast< Value >( 41 ); | |
196 | (*this)[8] = static_cast< Value >( 3 )/static_cast< Value >( 41 ); | |
197 | (*this)[9] = static_cast< Value >( 6 )/static_cast< Value >( 41 ); | |
198 | (*this)[10] = static_cast< Value >( 0 ); | |
199 | } | |
200 | }; | |
201 | ||
202 | template< class Value = double > | |
203 | struct rk78_coefficients_a12 : boost::array< Value , 12 > | |
204 | { | |
205 | rk78_coefficients_a12( void ) | |
206 | { | |
207 | (*this)[0] = static_cast< Value >( -1777 )/static_cast< Value >( 4100 ); | |
208 | (*this)[1] = static_cast< Value >( 0 ); | |
209 | (*this)[2] = static_cast< Value >( 0 ); | |
210 | (*this)[3] = static_cast< Value >( -341 )/static_cast< Value >( 164 ); | |
211 | (*this)[4] = static_cast< Value >( 4496 )/static_cast< Value >( 1025 ); | |
212 | (*this)[5] = static_cast< Value >( -289 )/static_cast< Value >( 82 ); | |
213 | (*this)[6] = static_cast< Value >( 2193 )/static_cast< Value >( 4100 ); | |
214 | (*this)[7] = static_cast< Value >( 51 )/static_cast< Value >( 82 ); | |
215 | (*this)[8] = static_cast< Value >( 33 )/static_cast< Value >( 164 ); | |
216 | (*this)[9] = static_cast< Value >( 12 )/static_cast< Value >( 41 ); | |
217 | (*this)[10] = static_cast< Value >( 0 ); | |
218 | (*this)[11] = static_cast< Value >( 1 ); | |
219 | } | |
220 | }; | |
221 | ||
222 | template< class Value = double > | |
223 | struct rk78_coefficients_b : boost::array< Value , 13 > | |
224 | { | |
225 | rk78_coefficients_b( void ) | |
226 | { | |
227 | (*this)[0] = static_cast< Value >( 0 ); | |
228 | (*this)[1] = static_cast< Value >( 0 ); | |
229 | (*this)[2] = static_cast< Value >( 0 ); | |
230 | (*this)[3] = static_cast< Value >( 0 ); | |
231 | (*this)[4] = static_cast< Value >( 0 ); | |
232 | (*this)[5] = static_cast< Value >( 34 )/static_cast<Value>( 105 ); | |
233 | (*this)[6] = static_cast< Value >( 9 )/static_cast<Value>( 35 ); | |
234 | (*this)[7] = static_cast< Value >( 9 )/static_cast<Value>( 35 ); | |
235 | (*this)[8] = static_cast< Value >( 9 )/static_cast<Value>( 280 ); | |
236 | (*this)[9] = static_cast< Value >( 9 )/static_cast<Value>( 280 ); | |
237 | (*this)[10] = static_cast< Value >( 0 ); | |
238 | (*this)[11] = static_cast< Value >( 41 )/static_cast<Value>( 840 ); | |
239 | (*this)[12] = static_cast< Value >( 41 )/static_cast<Value>( 840 ); | |
240 | } | |
241 | }; | |
242 | ||
243 | template< class Value = double > | |
244 | struct rk78_coefficients_db : boost::array< Value , 13 > | |
245 | { | |
246 | rk78_coefficients_db( void ) | |
247 | { | |
248 | (*this)[0] = static_cast< Value >( 0 ) - static_cast< Value >( 41 )/static_cast<Value>( 840 ); | |
249 | (*this)[1] = static_cast< Value >( 0 ); | |
250 | (*this)[2] = static_cast< Value >( 0 ); | |
251 | (*this)[3] = static_cast< Value >( 0 ); | |
252 | (*this)[4] = static_cast< Value >( 0 ); | |
253 | (*this)[5] = static_cast< Value >( 0 ); | |
254 | (*this)[6] = static_cast< Value >( 0 ); | |
255 | (*this)[7] = static_cast< Value >( 0 ); | |
256 | (*this)[8] = static_cast< Value >( 0 ); | |
257 | (*this)[9] = static_cast< Value >( 0 ); | |
258 | (*this)[10] = static_cast< Value >( 0 ) - static_cast< Value >( 41 )/static_cast<Value>( 840 ); | |
259 | (*this)[11] = static_cast< Value >( 41 )/static_cast<Value>( 840 ); | |
260 | (*this)[12] = static_cast< Value >( 41 )/static_cast<Value>( 840 ); | |
261 | } | |
262 | }; | |
263 | ||
264 | ||
265 | template< class Value = double > | |
266 | struct rk78_coefficients_c : boost::array< Value , 13 > | |
267 | { | |
268 | rk78_coefficients_c( void ) | |
269 | { | |
270 | (*this)[0] = static_cast< Value >( 0 ); | |
271 | (*this)[1] = static_cast< Value >( 2 )/static_cast< Value >( 27 ); | |
272 | (*this)[2] = static_cast< Value >( 1 )/static_cast< Value >( 9 ); | |
273 | (*this)[3] = static_cast< Value >( 1 )/static_cast<Value>( 6 ); | |
274 | (*this)[4] = static_cast< Value >( 5 )/static_cast<Value>( 12 ); | |
275 | (*this)[5] = static_cast< Value >( 1 )/static_cast<Value>( 2 ); | |
276 | (*this)[6] = static_cast< Value >( 5 )/static_cast<Value>( 6 ); | |
277 | (*this)[7] = static_cast< Value >( 1 )/static_cast<Value>( 6 ); | |
278 | (*this)[8] = static_cast< Value >( 2 )/static_cast<Value>( 3 ); | |
279 | (*this)[9] = static_cast< Value >( 1 )/static_cast<Value>( 3 ); | |
280 | (*this)[10] = static_cast< Value >( 1 ); | |
281 | (*this)[11] = static_cast< Value >( 0 ); | |
282 | (*this)[12] = static_cast< Value >( 1 ); | |
283 | } | |
284 | }; | |
285 | #endif // DOXYGEN_SKIP | |
286 | ||
287 | ||
288 | ||
289 | ||
290 | ||
291 | template< | |
292 | class State , | |
293 | class Value = double , | |
294 | class Deriv = State , | |
295 | class Time = Value , | |
296 | class Algebra = typename algebra_dispatcher< State >::algebra_type , | |
297 | class Operations = typename operations_dispatcher< State >::operations_type , | |
298 | class Resizer = initially_resizer | |
299 | > | |
300 | #ifndef DOXYGEN_SKIP | |
301 | class runge_kutta_fehlberg78 : public explicit_error_generic_rk< 13 , 8 , 8 , 7 , State , Value , Deriv , Time , | |
302 | Algebra , Operations , Resizer > | |
303 | #else | |
304 | class runge_kutta_fehlberg78 : public explicit_error_generic_rk | |
305 | #endif | |
306 | { | |
307 | ||
308 | public: | |
309 | #ifndef DOXYGEN_SKIP | |
310 | typedef explicit_error_generic_rk< 13 , 8 , 8 , 7 , State , Value , Deriv , Time , | |
311 | Algebra , Operations , Resizer > stepper_base_type; | |
312 | #endif | |
313 | typedef typename stepper_base_type::state_type state_type; | |
314 | typedef typename stepper_base_type::value_type value_type; | |
315 | typedef typename stepper_base_type::deriv_type deriv_type; | |
316 | typedef typename stepper_base_type::time_type time_type; | |
317 | typedef typename stepper_base_type::algebra_type algebra_type; | |
318 | typedef typename stepper_base_type::operations_type operations_type; | |
319 | typedef typename stepper_base_type::resizer_type resizer_type; | |
320 | ||
321 | #ifndef DOXYGEN_SKIP | |
322 | typedef typename stepper_base_type::stepper_type stepper_type; | |
323 | typedef typename stepper_base_type::wrapped_state_type wrapped_state_type; | |
324 | typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type; | |
325 | #endif // DOXYGEN_SKIP | |
326 | ||
327 | ||
328 | runge_kutta_fehlberg78( const algebra_type &algebra = algebra_type() ) : stepper_base_type( | |
329 | boost::fusion::make_vector( rk78_coefficients_a1<Value>() , rk78_coefficients_a2<Value>() , rk78_coefficients_a3<Value>() , | |
330 | rk78_coefficients_a4<Value>() , rk78_coefficients_a5<Value>() , rk78_coefficients_a6<Value>() , | |
331 | rk78_coefficients_a7<Value>() , rk78_coefficients_a8<Value>() , rk78_coefficients_a9<Value>() , | |
332 | rk78_coefficients_a10<Value>() , rk78_coefficients_a11<Value>() , rk78_coefficients_a12<Value>() ) , | |
333 | rk78_coefficients_b<Value>() , rk78_coefficients_db<Value>() , rk78_coefficients_c<Value>() , algebra ) | |
334 | { } | |
335 | }; | |
336 | ||
337 | ||
338 | ||
339 | /************* DOXYGEN *************/ | |
340 | ||
341 | /** | |
342 | * \class runge_kutta_fehlberg78 | |
343 | * \brief The Runge-Kutta Fehlberg 78 method. | |
344 | * | |
345 | * The Runge-Kutta Fehlberg 78 method is a standard method for high-precision applications. | |
346 | * The method is explicit and fulfills the Error Stepper concept. Step size control | |
347 | * is provided but continuous output is not available for this method. | |
348 | * | |
349 | * This class derives from explicit_error_stepper_base and inherits its interface via CRTP (current recurring template pattern). | |
350 | * Furthermore, it derivs from explicit_error_generic_rk which is a generic Runge-Kutta algorithm with error estimation. | |
351 | * For more details see explicit_error_stepper_base and explicit_error_generic_rk. | |
352 | * | |
353 | * \tparam State The state type. | |
354 | * \tparam Value The value type. | |
355 | * \tparam Deriv The type representing the time derivative of the state. | |
356 | * \tparam Time The time representing the independent variable - the time. | |
357 | * \tparam Algebra The algebra type. | |
358 | * \tparam Operations The operations type. | |
359 | * \tparam Resizer The resizer policy type. | |
360 | */ | |
361 | ||
362 | ||
363 | /** | |
364 | * \fn runge_kutta_fehlberg78::runge_kutta_fehlberg78( const algebra_type &algebra ) | |
365 | * \brief Constructs the runge_kutta_cash_fehlberg78 class. This constructor can be used as a default | |
366 | * constructor if the algebra has a default constructor. | |
367 | * \param algebra A copy of algebra is made and stored inside explicit_stepper_base. | |
368 | */ | |
369 | ||
370 | } | |
371 | } | |
372 | } | |
373 | ||
374 | #endif //BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_FEHLBERG87_HPP_INCLUDED |