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 |
375 | |