]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | /* |
2 | [auto_generated] | |
3 | boost/numeric/odeint/stepper/runge_kutta_cash_karp54_classic.hpp | |
4 | ||
5 | [begin_description] | |
6 | Classical implementation of the Runge-Kutta Cash-Karp 5(4) method. | |
7 | [end_description] | |
8 | ||
9 | Copyright 2010-2013 Mario Mulansky | |
10 | Copyright 2010-2013 Karsten Ahnert | |
11 | Copyright 2012 Christoph Koke | |
12 | ||
13 | Distributed under the Boost Software License, Version 1.0. | |
14 | (See accompanying file LICENSE_1_0.txt or | |
15 | copy at http://www.boost.org/LICENSE_1_0.txt) | |
16 | */ | |
17 | ||
18 | ||
19 | #ifndef BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_CASH_KARP54_CLASSIC_HPP_INCLUDED | |
20 | #define BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_CASH_KARP54_CLASSIC_HPP_INCLUDED | |
21 | ||
22 | ||
23 | #include <boost/numeric/odeint/util/bind.hpp> | |
24 | ||
25 | #include <boost/numeric/odeint/stepper/base/explicit_error_stepper_base.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 | #include <boost/numeric/odeint/stepper/stepper_categories.hpp> | |
31 | #include <boost/numeric/odeint/util/state_wrapper.hpp> | |
32 | #include <boost/numeric/odeint/util/is_resizeable.hpp> | |
33 | #include <boost/numeric/odeint/util/resizer.hpp> | |
34 | ||
35 | namespace boost { | |
36 | namespace numeric { | |
37 | namespace odeint { | |
38 | ||
39 | ||
40 | ||
41 | ||
42 | template< | |
43 | class State , | |
44 | class Value = double , | |
45 | class Deriv = State , | |
46 | class Time = Value , | |
47 | class Algebra = typename algebra_dispatcher< State >::algebra_type , | |
48 | class Operations = typename operations_dispatcher< State >::operations_type , | |
49 | class Resizer = initially_resizer | |
50 | > | |
51 | #ifndef DOXYGEN_SKIP | |
52 | class runge_kutta_cash_karp54_classic | |
53 | : public explicit_error_stepper_base< | |
54 | runge_kutta_cash_karp54_classic< State , Value , Deriv , Time , Algebra , Operations , Resizer > , | |
55 | 5 , 5 , 4 , State , Value , Deriv , Time , Algebra , Operations , Resizer > | |
56 | #else | |
57 | class runge_kutta_cash_karp54_classic : public explicit_error_stepper_base | |
58 | #endif | |
59 | { | |
60 | ||
61 | ||
62 | public : | |
63 | ||
64 | #ifndef DOXYGEN_SKIP | |
65 | typedef explicit_error_stepper_base< | |
66 | runge_kutta_cash_karp54_classic< State , Value , Deriv , Time , Algebra , Operations , Resizer > , | |
67 | 5 , 5 , 4 , State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_base_type; | |
68 | #else | |
69 | typedef explicit_error_stepper_base< runge_kutta_cash_karp54_classic< ... > , ... > stepper_base_type; | |
70 | #endif | |
71 | ||
72 | typedef typename stepper_base_type::state_type state_type; | |
73 | typedef typename stepper_base_type::value_type value_type; | |
74 | typedef typename stepper_base_type::deriv_type deriv_type; | |
75 | typedef typename stepper_base_type::time_type time_type; | |
76 | typedef typename stepper_base_type::algebra_type algebra_type; | |
77 | typedef typename stepper_base_type::operations_type operations_type; | |
78 | typedef typename stepper_base_type::resizer_type resizer_type; | |
79 | ||
80 | #ifndef DOXYGEN_SKIP | |
81 | typedef typename stepper_base_type::wrapped_state_type wrapped_state_type; | |
82 | typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type; | |
83 | typedef typename stepper_base_type::stepper_type stepper_type; | |
84 | #endif | |
85 | ||
86 | ||
87 | runge_kutta_cash_karp54_classic( const algebra_type &algebra = algebra_type() ) : stepper_base_type( algebra ) | |
88 | { } | |
89 | ||
90 | ||
91 | ||
92 | template< class System , class StateIn , class DerivIn , class StateOut , class Err > | |
93 | void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt , Err &xerr ) | |
94 | { | |
95 | const value_type c1 = static_cast<value_type> ( 37 ) / static_cast<value_type>( 378 ); | |
96 | const value_type c3 = static_cast<value_type> ( 250 ) / static_cast<value_type>( 621 ); | |
97 | const value_type c4 = static_cast<value_type> ( 125 ) / static_cast<value_type>( 594 ); | |
98 | const value_type c6 = static_cast<value_type> ( 512 ) / static_cast<value_type>( 1771 ); | |
99 | ||
100 | const value_type dc1 = c1 - static_cast<value_type> ( 2825 ) / static_cast<value_type>( 27648 ); | |
101 | const value_type dc3 = c3 - static_cast<value_type> ( 18575 ) / static_cast<value_type>( 48384 ); | |
102 | const value_type dc4 = c4 - static_cast<value_type> ( 13525 ) / static_cast<value_type>( 55296 ); | |
103 | const value_type dc5 = static_cast<value_type> ( -277 ) / static_cast<value_type>( 14336 ); | |
104 | const value_type dc6 = c6 - static_cast<value_type> ( 1 ) / static_cast<value_type> ( 4 ); | |
105 | ||
106 | do_step_impl( system , in , dxdt , t , out , dt ); | |
107 | ||
108 | //error estimate | |
109 | stepper_base_type::m_algebra.for_each6( xerr , dxdt , m_k3.m_v , m_k4.m_v , m_k5.m_v , m_k6.m_v , | |
110 | typename operations_type::template scale_sum5< time_type , time_type , time_type , time_type , time_type >( dt*dc1 , dt*dc3 , dt*dc4 , dt*dc5 , dt*dc6 )); | |
111 | ||
112 | } | |
113 | ||
114 | ||
115 | ||
116 | template< class System , class StateIn , class DerivIn , class StateOut > | |
117 | void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt ) | |
118 | { | |
119 | const value_type a2 = static_cast<value_type> ( 1 ) / static_cast<value_type> ( 5 ); | |
120 | const value_type a3 = static_cast<value_type> ( 3 ) / static_cast<value_type> ( 10 ); | |
121 | const value_type a4 = static_cast<value_type> ( 3 ) / static_cast<value_type> ( 5 ); | |
122 | const value_type a5 = static_cast<value_type> ( 1 ); | |
123 | const value_type a6 = static_cast<value_type> ( 7 ) / static_cast<value_type> ( 8 ); | |
124 | ||
125 | const value_type b21 = static_cast<value_type> ( 1 ) / static_cast<value_type> ( 5 ); | |
126 | const value_type b31 = static_cast<value_type> ( 3 ) / static_cast<value_type>( 40 ); | |
127 | const value_type b32 = static_cast<value_type> ( 9 ) / static_cast<value_type>( 40 ); | |
128 | const value_type b41 = static_cast<value_type> ( 3 ) / static_cast<value_type> ( 10 ); | |
129 | const value_type b42 = static_cast<value_type> ( -9 ) / static_cast<value_type> ( 10 ); | |
130 | const value_type b43 = static_cast<value_type> ( 6 ) / static_cast<value_type> ( 5 ); | |
131 | const value_type b51 = static_cast<value_type> ( -11 ) / static_cast<value_type>( 54 ); | |
132 | const value_type b52 = static_cast<value_type> ( 5 ) / static_cast<value_type> ( 2 ); | |
133 | const value_type b53 = static_cast<value_type> ( -70 ) / static_cast<value_type>( 27 ); | |
134 | const value_type b54 = static_cast<value_type> ( 35 ) / static_cast<value_type>( 27 ); | |
135 | const value_type b61 = static_cast<value_type> ( 1631 ) / static_cast<value_type>( 55296 ); | |
136 | const value_type b62 = static_cast<value_type> ( 175 ) / static_cast<value_type>( 512 ); | |
137 | const value_type b63 = static_cast<value_type> ( 575 ) / static_cast<value_type>( 13824 ); | |
138 | const value_type b64 = static_cast<value_type> ( 44275 ) / static_cast<value_type>( 110592 ); | |
139 | const value_type b65 = static_cast<value_type> ( 253 ) / static_cast<value_type>( 4096 ); | |
140 | ||
141 | const value_type c1 = static_cast<value_type> ( 37 ) / static_cast<value_type>( 378 ); | |
142 | const value_type c3 = static_cast<value_type> ( 250 ) / static_cast<value_type>( 621 ); | |
143 | const value_type c4 = static_cast<value_type> ( 125 ) / static_cast<value_type>( 594 ); | |
144 | const value_type c6 = static_cast<value_type> ( 512 ) / static_cast<value_type>( 1771 ); | |
145 | ||
146 | typename odeint::unwrap_reference< System >::type &sys = system; | |
147 | ||
148 | m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) ); | |
149 | ||
150 | //m_x1 = x + dt*b21*dxdt | |
151 | stepper_base_type::m_algebra.for_each3( m_x_tmp.m_v , in , dxdt , | |
152 | typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , dt*b21 ) ); | |
153 | ||
154 | sys( m_x_tmp.m_v , m_k2.m_v , t + dt*a2 ); | |
155 | // m_x_tmp = x + dt*b31*dxdt + dt*b32*m_x2 | |
156 | stepper_base_type::m_algebra.for_each4( m_x_tmp.m_v , in , dxdt , m_k2.m_v , | |
157 | typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt*b31 , dt*b32 )); | |
158 | ||
159 | sys( m_x_tmp.m_v , m_k3.m_v , t + dt*a3 ); | |
160 | // m_x_tmp = x + dt * (b41*dxdt + b42*m_x2 + b43*m_x3) | |
161 | stepper_base_type::m_algebra.for_each5( m_x_tmp.m_v , in , dxdt , m_k2.m_v , m_k3.m_v , | |
162 | typename operations_type::template scale_sum4< value_type , time_type , time_type , time_type >( 1.0 , dt*b41 , dt*b42 , dt*b43 )); | |
163 | ||
164 | sys( m_x_tmp.m_v, m_k4.m_v , t + dt*a4 ); | |
165 | stepper_base_type::m_algebra.for_each6( m_x_tmp.m_v , in , dxdt , m_k2.m_v , m_k3.m_v , m_k4.m_v , | |
166 | typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt*b51 , dt*b52 , dt*b53 , dt*b54 )); | |
167 | ||
168 | sys( m_x_tmp.m_v , m_k5.m_v , t + dt*a5 ); | |
169 | stepper_base_type::m_algebra.for_each7( m_x_tmp.m_v , in , dxdt , m_k2.m_v , m_k3.m_v , m_k4.m_v , m_k5.m_v , | |
170 | typename operations_type::template scale_sum6< value_type , time_type , time_type , time_type , time_type , time_type >( 1.0 , dt*b61 , dt*b62 , dt*b63 , dt*b64 , dt*b65 )); | |
171 | ||
172 | sys( m_x_tmp.m_v , m_k6.m_v , t + dt*a6 ); | |
173 | stepper_base_type::m_algebra.for_each6( out , in , dxdt , m_k3.m_v , m_k4.m_v , m_k6.m_v , | |
174 | typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt*c1 , dt*c3 , dt*c4 , dt*c6 )); | |
175 | ||
176 | } | |
177 | ||
178 | /** | |
179 | * \brief Adjust the size of all temporaries in the stepper manually. | |
180 | * \param x A state from which the size of the temporaries to be resized is deduced. | |
181 | */ | |
182 | template< class StateIn > | |
183 | void adjust_size( const StateIn &x ) | |
184 | { | |
185 | resize_impl( x ); | |
186 | stepper_base_type::adjust_size( x ); | |
187 | } | |
188 | ||
189 | private: | |
190 | ||
191 | template< class StateIn > | |
192 | bool resize_impl( const StateIn &x ) | |
193 | { | |
194 | bool resized = false; | |
195 | resized |= adjust_size_by_resizeability( m_x_tmp , x , typename is_resizeable<state_type>::type() ); | |
196 | resized |= adjust_size_by_resizeability( m_k2 , x , typename is_resizeable<deriv_type>::type() ); | |
197 | resized |= adjust_size_by_resizeability( m_k3 , x , typename is_resizeable<deriv_type>::type() ); | |
198 | resized |= adjust_size_by_resizeability( m_k4 , x , typename is_resizeable<deriv_type>::type() ); | |
199 | resized |= adjust_size_by_resizeability( m_k5 , x , typename is_resizeable<deriv_type>::type() ); | |
200 | resized |= adjust_size_by_resizeability( m_k6 , x , typename is_resizeable<deriv_type>::type() ); | |
201 | return resized; | |
202 | } | |
203 | ||
204 | ||
205 | wrapped_state_type m_x_tmp; | |
206 | wrapped_deriv_type m_k2, m_k3, m_k4, m_k5, m_k6; | |
207 | resizer_type m_resizer; | |
208 | ||
209 | }; | |
210 | ||
211 | ||
212 | ||
213 | /************ DOXYGEN *************/ | |
214 | ||
215 | /** | |
216 | * \class runge_kutta_cash_karp54_classic | |
217 | * \brief The Runge-Kutta Cash-Karp method implemented without the generic Runge-Kutta algorithm. | |
218 | * | |
219 | * The Runge-Kutta Cash-Karp method is one of the standard methods for | |
220 | * solving ordinary differential equations, see | |
221 | * <a href="http://en.wikipedia.org/wiki/Cash%E2%80%93Karp_method">en.wikipedia.org/wiki/Cash-Karp_method</a>. | |
222 | * The method is explicit and fulfills the Error Stepper concept. Step size control | |
223 | * is provided but continuous output is not available for this method. | |
224 | * | |
225 | * This class derives from explicit_error_stepper_base and inherits its interface via CRTP (current recurring | |
226 | * template pattern). This class implements the method directly, hence the generic Runge-Kutta algorithm is not used. | |
227 | * | |
228 | * \tparam State The state type. | |
229 | * \tparam Value The value type. | |
230 | * \tparam Deriv The type representing the time derivative of the state. | |
231 | * \tparam Time The time representing the independent variable - the time. | |
232 | * \tparam Algebra The algebra type. | |
233 | * \tparam Operations The operations type. | |
234 | * \tparam Resizer The resizer policy type. | |
235 | */ | |
236 | ||
237 | ||
238 | /** | |
239 | * \fn runge_kutta_cash_karp54_classic::runge_kutta_cash_karp54_classic( const algebra_type &algebra ) | |
240 | * \brief Constructs the runge_kutta_cash_karp54_classic class. This constructor can be used as a default | |
241 | * constructor if the algebra has a default constructor. | |
242 | * \param algebra A copy of algebra is made and stored inside explicit_stepper_base. | |
243 | */ | |
244 | ||
245 | ||
246 | /** | |
247 | * \fn runge_kutta_cash_karp54_classic::do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt , Err &xerr ) | |
248 | * \brief This method performs one step. The derivative `dxdt` of `in` at the time `t` is passed to the method. | |
249 | * | |
250 | * The result is updated out-of-place, hence the input is in `in` and the output in `out`. Futhermore, an | |
251 | * estimation of the error is stored in `xerr`. | |
252 | * Access to this step functionality is provided by explicit_error_stepper_base and | |
253 | * `do_step_impl` should not be called directly. | |
254 | ||
255 | * | |
256 | * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the | |
257 | * Simple System concept. | |
258 | * \param in The state of the ODE which should be solved. in is not modified in this method | |
259 | * \param dxdt The derivative of x at t. | |
260 | * \param t The value of the time, at which the step should be performed. | |
261 | * \param out The result of the step is written in out. | |
262 | * \param dt The step size. | |
263 | * \param xerr The result of the error estimation is written in xerr. | |
264 | */ | |
265 | ||
266 | /** | |
267 | * \fn runge_kutta_cash_karp54_classic::do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt ) | |
268 | * \brief This method performs one step. The derivative `dxdt` of `in` at the time `t` is passed to the method. | |
269 | * The result is updated out-of-place, hence the input is in `in` and the output in `out`. | |
270 | * Access to this step functionality is provided by explicit_error_stepper_base and | |
271 | * `do_step_impl` should not be called directly. | |
272 | * | |
273 | * \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the | |
274 | * Simple System concept. | |
275 | * \param in The state of the ODE which should be solved. in is not modified in this method | |
276 | * \param dxdt The derivative of x at t. | |
277 | * \param t The value of the time, at which the step should be performed. | |
278 | * \param out The result of the step is written in out. | |
279 | * \param dt The step size. | |
280 | */ | |
281 | ||
282 | } // odeint | |
283 | } // numeric | |
284 | } // boost | |
285 | ||
286 | ||
287 | ||
288 | ||
289 | #endif // BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_CASH_KARP54_CLASSIC_HPP_INCLUDED |