odeint/stepper/runge_kutta_fehlberg78.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/stepper/runge_kutta_fehlberg87.hpp
0004
0005  [begin_description]
0006  Implementation of the Runge-Kutta-Fehlberg stepper with the generic stepper.
0007  [end_description]
0008
0009  Copyright 2011-2013 Mario Mulansky
0010  Copyright 2012-2013 Karsten Ahnert
0011
0012  Distributed under the Boost Software License, Version 1.0.
0013  (See accompanying file LICENSE_1_0.txt or
0014  copy at http://www.boost.org/LICENSE_1_0.txt)
0015  */
0016
0017
0018 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_FEHLBERG87_HPP_INCLUDED
0019 #define BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_FEHLBERG87_HPP_INCLUDED
0020
0021
0022 #include <boost/fusion/container/vector.hpp>
0023 #include <boost/fusion/container/generation/make_vector.hpp>
0024
0025 #include <boost/numeric/odeint/stepper/explicit_error_generic_rk.hpp>
0026 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
0027 #include <boost/numeric/odeint/algebra/default_operations.hpp>
0028 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
0029 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
0030
0031 #include <boost/array.hpp>
0032
0033 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0034 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0035 #include <boost/numeric/odeint/util/resizer.hpp>
0036
0037
0038
0039
0040 namespace boost {
0041 namespace numeric {
0042 namespace odeint {
0043
0044
0045 #ifndef DOXYGEN_SKIP
0046 template< class Value = double >
0047 struct rk78_coefficients_a1 : boost::array< Value , 1 >
0048 {
0049     rk78_coefficients_a1( void )
0050             {
0051         (*this)[0] = static_cast< Value >( 2 )/static_cast< Value >( 27 );
0052             }
0053 };
0054
0055 template< class Value = double >
0056 struct rk78_coefficients_a2 : boost::array< Value , 2 >
0057 {
0058     rk78_coefficients_a2( void )
0059             {
0060         (*this)[0] = static_cast< Value >( 1 )/static_cast< Value >( 36 );
0061         (*this)[1] = static_cast< Value >( 1 )/static_cast< Value >( 12 );
0062             }
0063 };
0064
0065
0066 template< class Value = double >
0067 struct rk78_coefficients_a3 : boost::array< Value , 3 >
0068 {
0069     rk78_coefficients_a3( void )
0070             {
0071         (*this)[0] = static_cast< Value >( 1 )/static_cast< Value >( 24 );
0072         (*this)[1] = static_cast< Value >( 0 );
0073         (*this)[2] = static_cast< Value >( 1 )/static_cast< Value >( 8 );
0074             }
0075 };
0076
0077 template< class Value = double >
0078 struct rk78_coefficients_a4 : boost::array< Value , 4 >
0079 {
0080     rk78_coefficients_a4( void )
0081             {
0082         (*this)[0] = static_cast< Value >( 5 )/static_cast< Value >( 12 );
0083         (*this)[1] = static_cast< Value >( 0 );
0084         (*this)[2] = static_cast< Value >( -25 )/static_cast< Value >( 16 );
0085         (*this)[3] = static_cast< Value >( 25 )/static_cast< Value >( 16 );
0086             }
0087 };
0088
0089 template< class Value = double >
0090 struct rk78_coefficients_a5 : boost::array< Value , 5 >
0091 {
0092     rk78_coefficients_a5( void )
0093             {
0094         (*this)[0] = static_cast< Value >( 1 )/static_cast< Value >( 20 );
0095         (*this)[1] = static_cast< Value >( 0 );
0096         (*this)[2] = static_cast< Value >( 0 );
0097         (*this)[3] = static_cast< Value >( 1 )/static_cast< Value >( 4 );
0098         (*this)[4] = static_cast< Value >( 1 )/static_cast< Value >( 5 );
0099             }
0100 };
0101
0102
0103 template< class Value = double >
0104 struct rk78_coefficients_a6 : boost::array< Value , 6 >
0105 {
0106     rk78_coefficients_a6( void )
0107             {
0108         (*this)[0] = static_cast< Value >( -25 )/static_cast< Value >( 108 );
0109         (*this)[1] = static_cast< Value >( 0 );
0110         (*this)[2] = static_cast< Value >( 0 );
0111         (*this)[3] = static_cast< Value >( 125 )/static_cast< Value >( 108 );
0112         (*this)[4] = static_cast< Value >( -65 )/static_cast< Value >( 27 );
0113         (*this)[5] = static_cast< Value >( 125 )/static_cast< Value >( 54 );
0114             }
0115 };
0116
0117 template< class Value = double >
0118 struct rk78_coefficients_a7 : boost::array< Value , 7 >
0119 {
0120     rk78_coefficients_a7( void )
0121             {
0122         (*this)[0] = static_cast< Value >( 31 )/static_cast< Value >( 300 );
0123         (*this)[1] = static_cast< Value >( 0 );
0124         (*this)[2] = static_cast< Value >( 0 );
0125         (*this)[3] = static_cast< Value >( 0 );
0126         (*this)[4] = static_cast< Value >( 61 )/static_cast< Value >( 225 );
0127         (*this)[5] = static_cast< Value >( -2 )/static_cast< Value >( 9 );
0128         (*this)[6] = static_cast< Value >( 13 )/static_cast< Value >( 900 );
0129             }
0130 };
0131
0132 template< class Value = double >
0133 struct rk78_coefficients_a8 : boost::array< Value , 8 >
0134 {
0135     rk78_coefficients_a8( void )
0136             {
0137         (*this)[0] = static_cast< Value >( 2 );
0138         (*this)[1] = static_cast< Value >( 0 );
0139         (*this)[2] = static_cast< Value >( 0 );
0140         (*this)[3] = static_cast< Value >( -53 )/static_cast< Value >( 6 );
0141         (*this)[4] = static_cast< Value >( 704 )/static_cast< Value >( 45 );
0142         (*this)[5] = static_cast< Value >( -107 )/static_cast< Value >( 9 );
0143         (*this)[6] = static_cast< Value >( 67 )/static_cast< Value >( 90 );
0144         (*this)[7] = static_cast< Value >( 3 );
0145             }
0146 };
0147
0148 template< class Value = double >
0149 struct rk78_coefficients_a9 : boost::array< Value , 9 >
0150 {
0151     rk78_coefficients_a9( void )
0152             {
0153         (*this)[0] = static_cast< Value >( -91 )/static_cast< Value >( 108 );
0154         (*this)[1] = static_cast< Value >( 0 );
0155         (*this)[2] = static_cast< Value >( 0 );
0156         (*this)[3] = static_cast< Value >( 23 )/static_cast< Value >( 108 );
0157         (*this)[4] = static_cast< Value >( -976 )/static_cast< Value >( 135 );
0158         (*this)[5] = static_cast< Value >( 311 )/static_cast< Value >( 54 );
0159         (*this)[6] = static_cast< Value >( -19 )/static_cast< Value >( 60 );
0160         (*this)[7] = static_cast< Value >( 17 )/static_cast< Value >( 6 );
0161         (*this)[8] = static_cast< Value >( -1 )/static_cast< Value >( 12 );
0162             }
0163 };
0164
0165 template< class Value = double >
0166 struct rk78_coefficients_a10 : boost::array< Value , 10 >
0167 {
0168     rk78_coefficients_a10( void )
0169             {
0170         (*this)[0] = static_cast< Value >( 2383 )/static_cast< Value >( 4100 );
0171         (*this)[1] = static_cast< Value >( 0 );
0172         (*this)[2] = static_cast< Value >( 0 );
0173         (*this)[3] = static_cast< Value >( -341 )/static_cast< Value >( 164 );
0174         (*this)[4] = static_cast< Value >( 4496 )/static_cast< Value >( 1025 );
0175         (*this)[5] = static_cast< Value >( -301 )/static_cast< Value >( 82 );
0176         (*this)[6] = static_cast< Value >( 2133 )/static_cast< Value >( 4100 );
0177         (*this)[7] = static_cast< Value >( 45 )/static_cast< Value >( 82 );
0178         (*this)[8] = static_cast< Value >( 45 )/static_cast< Value >( 164 );
0179         (*this)[9] = static_cast< Value >( 18 )/static_cast< Value >( 41 );
0180             }
0181 };
0182
0183 template< class Value = double >
0184 struct rk78_coefficients_a11 : boost::array< Value , 11 >
0185 {
0186     rk78_coefficients_a11( void )
0187             {
0188         (*this)[0] = static_cast< Value >( 3 )/static_cast< Value >( 205 );
0189         (*this)[1] = static_cast< Value >( 0 );
0190         (*this)[2] = static_cast< Value >( 0 );
0191         (*this)[3] = static_cast< Value >( 0 );
0192         (*this)[4] = static_cast< Value >( 0 );
0193         (*this)[5] = static_cast< Value >( -6 )/static_cast< Value >( 41 );
0194         (*this)[6] = static_cast< Value >( -3 )/static_cast< Value >( 205 );
0195         (*this)[7] = static_cast< Value >( -3 )/static_cast< Value >( 41 );
0196         (*this)[8] = static_cast< Value >( 3 )/static_cast< Value >( 41 );
0197         (*this)[9] = static_cast< Value >( 6 )/static_cast< Value >( 41 );
0198         (*this)[10] = static_cast< Value >( 0 );
0199             }
0200 };
0201
0202 template< class Value = double >
0203 struct rk78_coefficients_a12 : boost::array< Value , 12 >
0204 {
0205     rk78_coefficients_a12( void )
0206             {
0207         (*this)[0] = static_cast< Value >( -1777 )/static_cast< Value >( 4100 );
0208         (*this)[1] = static_cast< Value >( 0 );
0209         (*this)[2] = static_cast< Value >( 0 );
0210         (*this)[3] = static_cast< Value >( -341 )/static_cast< Value >( 164 );
0211         (*this)[4] = static_cast< Value >( 4496 )/static_cast< Value >( 1025 );
0212         (*this)[5] = static_cast< Value >( -289 )/static_cast< Value >( 82 );
0213         (*this)[6] = static_cast< Value >( 2193 )/static_cast< Value >( 4100 );
0214         (*this)[7] = static_cast< Value >( 51 )/static_cast< Value >( 82 );
0215         (*this)[8] = static_cast< Value >( 33 )/static_cast< Value >( 164 );
0216         (*this)[9] = static_cast< Value >( 12 )/static_cast< Value >( 41 );
0217         (*this)[10] = static_cast< Value >( 0 );
0218         (*this)[11] = static_cast< Value >( 1 );
0219             }
0220 };
0221
0222 template< class Value = double >
0223 struct rk78_coefficients_b : boost::array< Value , 13 >
0224 {
0225     rk78_coefficients_b( void )
0226             {
0227         (*this)[0] = static_cast< Value >( 0 );
0228         (*this)[1] = static_cast< Value >( 0 );
0229         (*this)[2] = static_cast< Value >( 0 );
0230         (*this)[3] = static_cast< Value >( 0 );
0231         (*this)[4] = static_cast< Value >( 0 );
0232         (*this)[5] = static_cast< Value >( 34 )/static_cast<Value>( 105 );
0233         (*this)[6] = static_cast< Value >( 9 )/static_cast<Value>( 35 );
0234         (*this)[7] = static_cast< Value >( 9 )/static_cast<Value>( 35 );
0235         (*this)[8] = static_cast< Value >( 9 )/static_cast<Value>( 280 );
0236         (*this)[9] = static_cast< Value >( 9 )/static_cast<Value>( 280 );
0237         (*this)[10] = static_cast< Value >( 0 );
0238         (*this)[11] = static_cast< Value >( 41 )/static_cast<Value>( 840 );
0239         (*this)[12] = static_cast< Value >( 41 )/static_cast<Value>( 840 );
0240             }
0241 };
0242
0243 template< class Value = double >
0244 struct rk78_coefficients_db : boost::array< Value , 13 >
0245 {
0246     rk78_coefficients_db( void )
0247             {
0248         (*this)[0] = static_cast< Value >( 0 ) - static_cast< Value >( 41 )/static_cast<Value>( 840 );
0249         (*this)[1] = static_cast< Value >( 0 );
0250         (*this)[2] = static_cast< Value >( 0 );
0251         (*this)[3] = static_cast< Value >( 0 );
0252         (*this)[4] = static_cast< Value >( 0 );
0253         (*this)[5] = static_cast< Value >( 0 );
0254         (*this)[6] = static_cast< Value >( 0 );
0255         (*this)[7] = static_cast< Value >( 0 );
0256         (*this)[8] = static_cast< Value >( 0 );
0257         (*this)[9] = static_cast< Value >( 0 );
0258         (*this)[10] = static_cast< Value >( 0 ) - static_cast< Value >( 41 )/static_cast<Value>( 840 );
0259         (*this)[11] = static_cast< Value >( 41 )/static_cast<Value>( 840 );
0260         (*this)[12] = static_cast< Value >( 41 )/static_cast<Value>( 840 );
0261             }
0262 };
0263
0264
0265 template< class Value = double >
0266 struct rk78_coefficients_c : boost::array< Value , 13 >
0267 {
0268     rk78_coefficients_c( void )
0269             {
0270         (*this)[0] = static_cast< Value >( 0 );
0271         (*this)[1] = static_cast< Value >( 2 )/static_cast< Value >( 27 );
0272         (*this)[2] = static_cast< Value >( 1 )/static_cast< Value >( 9 );
0273         (*this)[3] = static_cast< Value >( 1 )/static_cast<Value>( 6 );
0274         (*this)[4] = static_cast< Value >( 5 )/static_cast<Value>( 12 );
0275         (*this)[5] = static_cast< Value >( 1 )/static_cast<Value>( 2 );
0276         (*this)[6] = static_cast< Value >( 5 )/static_cast<Value>( 6 );
0277         (*this)[7] = static_cast< Value >( 1 )/static_cast<Value>( 6 );
0278         (*this)[8] = static_cast< Value >( 2 )/static_cast<Value>( 3 );
0279         (*this)[9] = static_cast< Value >( 1 )/static_cast<Value>( 3 );
0280         (*this)[10] = static_cast< Value >( 1 );
0281         (*this)[11] = static_cast< Value >( 0 );
0282         (*this)[12] = static_cast< Value >( 1 );
0283             }
0284 };
0285 #endif // DOXYGEN_SKIP
0286
0287
0288
0289
0290
0291 template<
0292 class State ,
0293 class Value = double ,
0294 class Deriv = State ,
0295 class Time = Value ,
0296 class Algebra = typename algebra_dispatcher< State >::algebra_type ,
0297 class Operations = typename operations_dispatcher< State >::operations_type ,
0298 class Resizer = initially_resizer
0299 >
0300 #ifndef DOXYGEN_SKIP
0301 class runge_kutta_fehlberg78 : public explicit_error_generic_rk< 13 , 8 , 8 , 7 , State , Value , Deriv , Time ,
0302 Algebra , Operations , Resizer >
0303 #else
0304 class runge_kutta_fehlberg78 : public explicit_error_generic_rk
0305 #endif
0306 {
0307
0308 public:
0309 #ifndef DOXYGEN_SKIP
0310     typedef explicit_error_generic_rk< 13 , 8 , 8 , 7 , State , Value , Deriv , Time ,
0311             Algebra , Operations , Resizer > stepper_base_type;
0312 #endif
0313     typedef typename stepper_base_type::state_type state_type;
0314     typedef typename stepper_base_type::value_type value_type;
0315     typedef typename stepper_base_type::deriv_type deriv_type;
0316     typedef typename stepper_base_type::time_type time_type;
0317     typedef typename stepper_base_type::algebra_type algebra_type;
0318     typedef typename stepper_base_type::operations_type operations_type;
0319     typedef typename stepper_base_type::resizer_type resizer_type;
0320
0321     #ifndef DOXYGEN_SKIP
0322     typedef typename stepper_base_type::stepper_type stepper_type;
0323     typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
0324     typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
0325     #endif // DOXYGEN_SKIP
0326
0327
0328     runge_kutta_fehlberg78( const algebra_type &algebra = algebra_type() ) : stepper_base_type(
0329             boost::fusion::make_vector( rk78_coefficients_a1<Value>() , rk78_coefficients_a2<Value>() , rk78_coefficients_a3<Value>() ,
0330                     rk78_coefficients_a4<Value>() , rk78_coefficients_a5<Value>() , rk78_coefficients_a6<Value>() ,
0331                     rk78_coefficients_a7<Value>() , rk78_coefficients_a8<Value>() , rk78_coefficients_a9<Value>() ,
0332                     rk78_coefficients_a10<Value>() , rk78_coefficients_a11<Value>() , rk78_coefficients_a12<Value>() ) ,
0333             rk78_coefficients_b<Value>() , rk78_coefficients_db<Value>() , rk78_coefficients_c<Value>() , algebra )
0334     { }
0335 };
0336
0337
0338
0339 /************* DOXYGEN *************/
0340
0341 /**
0342  * \class runge_kutta_fehlberg78
0343  * \brief The Runge-Kutta Fehlberg 78 method.
0344  *
0345  * The Runge-Kutta Fehlberg 78 method is a standard method for high-precision applications.
0346  * The method is explicit and fulfills the Error Stepper concept. Step size control
0347  * is provided but continuous output is not available for this method.
0348  *
0349  * This class derives from explicit_error_stepper_base and inherits its interface via CRTP (current recurring template pattern).
0350  * Furthermore, it derivs from explicit_error_generic_rk which is a generic Runge-Kutta algorithm with error estimation.
0351  * For more details see explicit_error_stepper_base and explicit_error_generic_rk.
0352  *
0353  * \tparam State The state type.
0354  * \tparam Value The value type.
0355  * \tparam Deriv The type representing the time derivative of the state.
0356  * \tparam Time The time representing the independent variable - the time.
0357  * \tparam Algebra The algebra type.
0358  * \tparam Operations The operations type.
0359  * \tparam Resizer The resizer policy type.
0360  */
0361
0362
0363     /**
0364      * \fn runge_kutta_fehlberg78::runge_kutta_fehlberg78( const algebra_type &algebra )
0365      * \brief Constructs the runge_kutta_cash_fehlberg78 class. This constructor can be used as a default
0366      * constructor if the algebra has a default constructor.
0367      * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
0368      */
0369
0370 }
0371 }
0372 }
0373
0374 #endif //BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_FEHLBERG87_HPP_INCLUDED