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 { rk78_coefficients_a1boost::numeric::odeint::rk78_coefficients_a149 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 { rk78_coefficients_a2boost::numeric::odeint::rk78_coefficients_a258 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 { rk78_coefficients_a3boost::numeric::odeint::rk78_coefficients_a369 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 { rk78_coefficients_a4boost::numeric::odeint::rk78_coefficients_a480 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 { rk78_coefficients_a5boost::numeric::odeint::rk78_coefficients_a592 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 { rk78_coefficients_a6boost::numeric::odeint::rk78_coefficients_a6106 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 { rk78_coefficients_a7boost::numeric::odeint::rk78_coefficients_a7120 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 { rk78_coefficients_a8boost::numeric::odeint::rk78_coefficients_a8135 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 { rk78_coefficients_a9boost::numeric::odeint::rk78_coefficients_a9151 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 { rk78_coefficients_a10boost::numeric::odeint::rk78_coefficients_a10168 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 { rk78_coefficients_a11boost::numeric::odeint::rk78_coefficients_a11186 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 { rk78_coefficients_a12boost::numeric::odeint::rk78_coefficients_a12205 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 { rk78_coefficients_bboost::numeric::odeint::rk78_coefficients_b225 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 { rk78_coefficients_dbboost::numeric::odeint::rk78_coefficients_db246 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 { rk78_coefficients_cboost::numeric::odeint::rk78_coefficients_c268 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 runge_kutta_fehlberg78(const algebra_type & algebra=algebra_type ())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