runge_kutta_cash_karp54.hpp 8.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231
  1. /*
  2. [auto_generated]
  3. boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp
  4. [begin_description]
  5. Implementation of the Runge Kutta Cash Karp 5(4) method. It uses the generic error stepper.
  6. [end_description]
  7. Copyright 2011-2013 Mario Mulansky
  8. Copyright 2011-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_CASH_KARP54_HPP_INCLUDED
  14. #define BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_CASH_KARP54_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/numeric/odeint/util/state_wrapper.hpp>
  23. #include <boost/numeric/odeint/util/is_resizeable.hpp>
  24. #include <boost/numeric/odeint/util/resizer.hpp>
  25. #include <boost/array.hpp>
  26. namespace boost {
  27. namespace numeric {
  28. namespace odeint {
  29. #ifndef DOXYGEN_SKIP
  30. template< class Value = double >
  31. struct rk54_ck_coefficients_a1 : boost::array< Value , 1 >
  32. {
  33. rk54_ck_coefficients_a1( void )
  34. {
  35. (*this)[0] = static_cast< Value >( 1 )/static_cast< Value >( 5 );
  36. }
  37. };
  38. template< class Value = double >
  39. struct rk54_ck_coefficients_a2 : boost::array< Value , 2 >
  40. {
  41. rk54_ck_coefficients_a2( void )
  42. {
  43. (*this)[0] = static_cast<Value>( 3 )/static_cast<Value>( 40 );
  44. (*this)[1] = static_cast<Value>( 9 )/static_cast<Value>( 40 );
  45. }
  46. };
  47. template< class Value = double >
  48. struct rk54_ck_coefficients_a3 : boost::array< Value , 3 >
  49. {
  50. rk54_ck_coefficients_a3( void )
  51. {
  52. (*this)[0] = static_cast<Value>( 3 )/static_cast<Value>( 10 );
  53. (*this)[1] = static_cast<Value>( -9 )/static_cast<Value>( 10 );
  54. (*this)[2] = static_cast<Value>( 6 )/static_cast<Value>( 5 );
  55. }
  56. };
  57. template< class Value = double >
  58. struct rk54_ck_coefficients_a4 : boost::array< Value , 4 >
  59. {
  60. rk54_ck_coefficients_a4( void )
  61. {
  62. (*this)[0] = static_cast<Value>( -11 )/static_cast<Value>( 54 );
  63. (*this)[1] = static_cast<Value>( 5 )/static_cast<Value>( 2 );
  64. (*this)[2] = static_cast<Value>( -70 )/static_cast<Value>( 27 );
  65. (*this)[3] = static_cast<Value>( 35 )/static_cast<Value>( 27 );
  66. }
  67. };
  68. template< class Value = double >
  69. struct rk54_ck_coefficients_a5 : boost::array< Value , 5 >
  70. {
  71. rk54_ck_coefficients_a5( void )
  72. {
  73. (*this)[0] = static_cast<Value>( 1631 )/static_cast<Value>( 55296 );
  74. (*this)[1] = static_cast<Value>( 175 )/static_cast<Value>( 512 );
  75. (*this)[2] = static_cast<Value>( 575 )/static_cast<Value>( 13824 );
  76. (*this)[3] = static_cast<Value>( 44275 )/static_cast<Value>( 110592 );
  77. (*this)[4] = static_cast<Value>( 253 )/static_cast<Value>( 4096 );
  78. }
  79. };
  80. template< class Value = double >
  81. struct rk54_ck_coefficients_b : boost::array< Value , 6 >
  82. {
  83. rk54_ck_coefficients_b( void )
  84. {
  85. (*this)[0] = static_cast<Value>( 37 )/static_cast<Value>( 378 );
  86. (*this)[1] = static_cast<Value>( 0 );
  87. (*this)[2] = static_cast<Value>( 250 )/static_cast<Value>( 621 );
  88. (*this)[3] = static_cast<Value>( 125 )/static_cast<Value>( 594 );
  89. (*this)[4] = static_cast<Value>( 0 );
  90. (*this)[5] = static_cast<Value>( 512 )/static_cast<Value>( 1771 );
  91. }
  92. };
  93. template< class Value = double >
  94. struct rk54_ck_coefficients_db : boost::array< Value , 6 >
  95. {
  96. rk54_ck_coefficients_db( void )
  97. {
  98. (*this)[0] = static_cast<Value>( 37 )/static_cast<Value>( 378 ) - static_cast<Value>( 2825 )/static_cast<Value>( 27648 );
  99. (*this)[1] = static_cast<Value>( 0 );
  100. (*this)[2] = static_cast<Value>( 250 )/static_cast<Value>( 621 ) - static_cast<Value>( 18575 )/static_cast<Value>( 48384 );
  101. (*this)[3] = static_cast<Value>( 125 )/static_cast<Value>( 594 ) - static_cast<Value>( 13525 )/static_cast<Value>( 55296 );
  102. (*this)[4] = static_cast<Value>( -277 )/static_cast<Value>( 14336 );
  103. (*this)[5] = static_cast<Value>( 512 )/static_cast<Value>( 1771 ) - static_cast<Value>( 1 )/static_cast<Value>( 4 );
  104. }
  105. };
  106. template< class Value = double >
  107. struct rk54_ck_coefficients_c : boost::array< Value , 6 >
  108. {
  109. rk54_ck_coefficients_c( void )
  110. {
  111. (*this)[0] = static_cast<Value>(0);
  112. (*this)[1] = static_cast<Value>( 1 )/static_cast<Value>( 5 );
  113. (*this)[2] = static_cast<Value>( 3 )/static_cast<Value>( 10 );
  114. (*this)[3] = static_cast<Value>( 3 )/static_cast<Value>( 5 );
  115. (*this)[4] = static_cast<Value>( 1 );
  116. (*this)[5] = static_cast<Value>( 7 )/static_cast<Value>( 8 );
  117. }
  118. };
  119. #endif
  120. template<
  121. class State ,
  122. class Value = double ,
  123. class Deriv = State ,
  124. class Time = Value ,
  125. class Algebra = typename algebra_dispatcher< State >::algebra_type ,
  126. class Operations = typename operations_dispatcher< State >::operations_type ,
  127. class Resizer = initially_resizer
  128. >
  129. #ifndef DOXYGEN_SKIP
  130. class runge_kutta_cash_karp54 : public explicit_error_generic_rk< 6 , 5 , 5 , 4 ,
  131. State , Value , Deriv , Time , Algebra , Operations , Resizer >
  132. #else
  133. class runge_kutta_cash_karp54 : public explicit_error_generic_rk
  134. #endif
  135. {
  136. public:
  137. #ifndef DOXYGEN_SKIP
  138. typedef explicit_error_generic_rk< 6 , 5 , 5 , 4 , State , Value , Deriv , Time ,
  139. Algebra , Operations , Resizer > stepper_base_type;
  140. #endif
  141. typedef typename stepper_base_type::state_type state_type;
  142. typedef typename stepper_base_type::value_type value_type;
  143. typedef typename stepper_base_type::deriv_type deriv_type;
  144. typedef typename stepper_base_type::time_type time_type;
  145. typedef typename stepper_base_type::algebra_type algebra_type;
  146. typedef typename stepper_base_type::operations_type operations_type;
  147. typedef typename stepper_base_type::resizer_type resizer_typ;
  148. #ifndef DOXYGEN_SKIP
  149. typedef typename stepper_base_type::stepper_type stepper_type;
  150. typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
  151. typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
  152. #endif
  153. runge_kutta_cash_karp54( const algebra_type &algebra = algebra_type() ) : stepper_base_type(
  154. boost::fusion::make_vector( rk54_ck_coefficients_a1<Value>() ,
  155. rk54_ck_coefficients_a2<Value>() ,
  156. rk54_ck_coefficients_a3<Value>() ,
  157. rk54_ck_coefficients_a4<Value>() ,
  158. rk54_ck_coefficients_a5<Value>() ) ,
  159. rk54_ck_coefficients_b<Value>() , rk54_ck_coefficients_db<Value>() , rk54_ck_coefficients_c<Value>() ,
  160. algebra )
  161. { }
  162. };
  163. /********** DOXYGEN **********/
  164. /**
  165. * \class runge_kutta_cash_karp54
  166. * \brief The Runge-Kutta Cash-Karp method.
  167. *
  168. * The Runge-Kutta Cash-Karp method is one of the standard methods for
  169. * solving ordinary differential equations, see
  170. * <a href="http://en.wikipedia.org/wiki/Cash%E2%80%93Karp_methods">en.wikipedia.org/wiki/Cash-Karp_methods</a>.
  171. * The method is explicit and fulfills the Error Stepper concept. Step size control
  172. * is provided but continuous output is not available for this method.
  173. *
  174. * This class derives from explicit_error_stepper_base and inherits its interface via CRTP (current recurring template pattern).
  175. * Furthermore, it derivs from explicit_error_generic_rk which is a generic Runge-Kutta algorithm with error estimation.
  176. * For more details see explicit_error_stepper_base and explicit_error_generic_rk.
  177. *
  178. * \tparam State The state type.
  179. * \tparam Value The value type.
  180. * \tparam Deriv The type representing the time derivative of the state.
  181. * \tparam Time The time representing the independent variable - the time.
  182. * \tparam Algebra The algebra type.
  183. * \tparam Operations The operations type.
  184. * \tparam Resizer The resizer policy type.
  185. */
  186. /**
  187. * \fn runge_kutta_cash_karp54::runge_kutta_cash_karp54( const algebra_type &algebra )
  188. * \brief Constructs the runge_kutta_cash_karp54 class. This constructor can be used as a default
  189. * constructor if the algebra has a default constructor.
  190. * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
  191. */
  192. }
  193. }
  194. }
  195. #endif // BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_CASH_KARP54_HPP_INCLUDED