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
40namespace boost {
41namespace numeric {
42namespace odeint {
43
44
45#ifndef DOXYGEN_SKIP
46template< class Value = double >
47struct 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
55template< class Value = double >
56struct 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
66template< class Value = double >
67struct 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
77template< class Value = double >
78struct 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
89template< class Value = double >
90struct 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
103template< class Value = double >
104struct 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
117template< class Value = double >
118struct 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
132template< class Value = double >
133struct 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
148template< class Value = double >
149struct 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
165template< class Value = double >
166struct 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
183template< class Value = double >
184struct 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
202template< class Value = double >
203struct 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
222template< class Value = double >
223struct 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
243template< class Value = double >
244struct 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
265template< class Value = double >
266struct 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
291template<
292class State ,
293class Value = double ,
294class Deriv = State ,
295class Time = Value ,
296class Algebra = typename algebra_dispatcher< State >::algebra_type ,
297class Operations = typename operations_dispatcher< State >::operations_type ,
298class Resizer = initially_resizer
299>
300#ifndef DOXYGEN_SKIP
301class runge_kutta_fehlberg78 : public explicit_error_generic_rk< 13 , 8 , 8 , 7 , State , Value , Deriv , Time ,
302Algebra , Operations , Resizer >
303#else
304class runge_kutta_fehlberg78 : public explicit_error_generic_rk
305#endif
306{
307
308public:
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

source code of boost/boost/numeric/odeint/stepper/runge_kutta_fehlberg78.hpp