runge_kutta_fehlberg78.hpp 14 KB

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