runge_kutta_fehlberg78.hpp 15 KB

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