]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | /* |
2 | [auto_generated] | |
3 | boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp | |
4 | ||
5 | [begin_description] | |
6 | Implementation of the Dense-output stepper for all steppers. Note, that this class does | |
7 | not computes the result but serves as an interface. | |
8 | [end_description] | |
9 | ||
10 | Copyright 2011-2013 Karsten Ahnert | |
11 | Copyright 2011-2015 Mario Mulansky | |
12 | Copyright 2012 Christoph Koke | |
13 | ||
14 | Distributed under the Boost Software License, Version 1.0. | |
15 | (See accompanying file LICENSE_1_0.txt or | |
16 | copy at http://www.boost.org/LICENSE_1_0.txt) | |
17 | */ | |
18 | ||
19 | ||
20 | #ifndef BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED | |
21 | #define BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED | |
22 | ||
23 | ||
24 | #include <utility> | |
25 | #include <stdexcept> | |
26 | ||
27 | #include <boost/throw_exception.hpp> | |
28 | ||
29 | #include <boost/numeric/odeint/util/bind.hpp> | |
30 | ||
31 | #include <boost/numeric/odeint/util/copy.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 | #include <boost/numeric/odeint/stepper/controlled_step_result.hpp> | |
38 | #include <boost/numeric/odeint/stepper/stepper_categories.hpp> | |
39 | ||
40 | #include <boost/numeric/odeint/integrate/max_step_checker.hpp> | |
41 | ||
42 | namespace boost { | |
43 | namespace numeric { | |
44 | namespace odeint { | |
45 | ||
46 | template< class Stepper , class StepperCategory = typename Stepper::stepper_category > | |
47 | class dense_output_runge_kutta; | |
48 | ||
49 | ||
50 | /** | |
51 | * \brief The class representing dense-output Runge-Kutta steppers. | |
52 | * \note In this stepper, the initialize method has to be called before using | |
53 | * the do_step method. | |
54 | * | |
55 | * The dense-output functionality allows to interpolate the solution between | |
56 | * subsequent integration points using intermediate results obtained during the | |
57 | * computation. This version works based on a normal stepper without step-size | |
58 | * control. | |
59 | * | |
60 | * | |
61 | * \tparam Stepper The stepper type of the underlying algorithm. | |
62 | */ | |
63 | template< class Stepper > | |
64 | class dense_output_runge_kutta< Stepper , stepper_tag > | |
65 | { | |
66 | ||
67 | public: | |
68 | ||
69 | /* | |
70 | * We do not need all typedefs. | |
71 | */ | |
72 | typedef Stepper stepper_type; | |
73 | typedef typename stepper_type::state_type state_type; | |
74 | typedef typename stepper_type::wrapped_state_type wrapped_state_type; | |
75 | typedef typename stepper_type::value_type value_type; | |
76 | typedef typename stepper_type::deriv_type deriv_type; | |
77 | typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type; | |
78 | typedef typename stepper_type::time_type time_type; | |
79 | typedef typename stepper_type::algebra_type algebra_type; | |
80 | typedef typename stepper_type::operations_type operations_type; | |
81 | typedef typename stepper_type::resizer_type resizer_type; | |
82 | typedef dense_output_stepper_tag stepper_category; | |
83 | typedef dense_output_runge_kutta< Stepper > dense_output_stepper_type; | |
84 | ||
85 | ||
86 | /** | |
87 | * \brief Constructs the dense_output_runge_kutta class. An instance of the | |
88 | * underlying stepper can be provided. | |
89 | * \param stepper An instance of the underlying stepper. | |
90 | */ | |
91 | dense_output_runge_kutta( const stepper_type &stepper = stepper_type() ) | |
92 | : m_stepper( stepper ) , m_resizer() , | |
93 | m_x1() , m_x2() , m_current_state_x1( true ) , | |
94 | m_t() , m_t_old() , m_dt() | |
95 | { } | |
96 | ||
97 | ||
98 | /** | |
99 | * \brief Initializes the stepper. Has to be called before do_step can be | |
100 | * used to set the initial conditions and the step size. | |
101 | * \param x0 The initial state of the ODE which should be solved. | |
102 | * \param t0 The initial time, at which the step should be performed. | |
103 | * \param dt0 The step size. | |
104 | */ | |
105 | template< class StateType > | |
106 | void initialize( const StateType &x0 , time_type t0 , time_type dt0 ) | |
107 | { | |
108 | m_resizer.adjust_size( x0 , detail::bind( &dense_output_stepper_type::template resize_impl< StateType > , detail::ref( *this ) , detail::_1 ) ); | |
109 | boost::numeric::odeint::copy( x0 , get_current_state() ); | |
110 | m_t = t0; | |
111 | m_dt = dt0; | |
112 | } | |
113 | ||
114 | /** | |
115 | * \brief Does one time step. | |
116 | * \note initialize has to be called before using this method to set the | |
117 | * initial conditions x,t and the stepsize. | |
118 | * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the | |
119 | * Simple System concept. | |
120 | * \return Pair with start and end time of the integration step. | |
121 | */ | |
122 | template< class System > | |
123 | std::pair< time_type , time_type > do_step( System system ) | |
124 | { | |
125 | m_stepper.do_step( system , get_current_state() , m_t , get_old_state() , m_dt ); | |
126 | m_t_old = m_t; | |
127 | m_t += m_dt; | |
128 | toggle_current_state(); | |
92f5a8d4 | 129 | return std::make_pair( m_t_old , m_t ); |
7c673cae FG |
130 | } |
131 | ||
132 | /* | |
133 | * The next two overloads are needed to solve the forwarding problem | |
134 | */ | |
135 | ||
136 | /** | |
137 | * \brief Calculates the solution at an intermediate point. | |
138 | * \param t The time at which the solution should be calculated, has to be | |
139 | * in the current time interval. | |
140 | * \param x The output variable where the result is written into. | |
141 | */ | |
142 | template< class StateOut > | |
143 | void calc_state( time_type t , StateOut &x ) const | |
144 | { | |
145 | if( t == current_time() ) | |
146 | { | |
147 | boost::numeric::odeint::copy( get_current_state() , x ); | |
148 | } | |
149 | m_stepper.calc_state( x , t , get_old_state() , m_t_old , get_current_state() , m_t ); | |
150 | } | |
151 | ||
152 | /** | |
153 | * \brief Calculates the solution at an intermediate point. Solves the forwarding problem | |
154 | * \param t The time at which the solution should be calculated, has to be | |
155 | * in the current time interval. | |
156 | * \param x The output variable where the result is written into, can be a boost range. | |
157 | */ | |
158 | template< class StateOut > | |
159 | void calc_state( time_type t , const StateOut &x ) const | |
160 | { | |
161 | m_stepper.calc_state( x , t , get_old_state() , m_t_old , get_current_state() , m_t ); | |
162 | } | |
163 | ||
164 | /** | |
165 | * \brief Adjust the size of all temporaries in the stepper manually. | |
166 | * \param x A state from which the size of the temporaries to be resized is deduced. | |
167 | */ | |
168 | template< class StateType > | |
169 | void adjust_size( const StateType &x ) | |
170 | { | |
171 | resize_impl( x ); | |
172 | m_stepper.stepper().resize( x ); | |
173 | } | |
174 | ||
175 | /** | |
176 | * \brief Returns the current state of the solution. | |
177 | * \return The current state of the solution x(t). | |
178 | */ | |
179 | const state_type& current_state( void ) const | |
180 | { | |
181 | return get_current_state(); | |
182 | } | |
183 | ||
184 | /** | |
185 | * \brief Returns the current time of the solution. | |
186 | * \return The current time of the solution t. | |
187 | */ | |
188 | time_type current_time( void ) const | |
189 | { | |
190 | return m_t; | |
191 | } | |
192 | ||
193 | /** | |
194 | * \brief Returns the last state of the solution. | |
195 | * \return The last state of the solution x(t-dt). | |
196 | */ | |
197 | const state_type& previous_state( void ) const | |
198 | { | |
199 | return get_old_state(); | |
200 | } | |
201 | ||
202 | /** | |
203 | * \brief Returns the last time of the solution. | |
204 | * \return The last time of the solution t-dt. | |
205 | */ | |
206 | time_type previous_time( void ) const | |
207 | { | |
208 | return m_t_old; | |
209 | } | |
210 | ||
211 | /** | |
212 | * \brief Returns the current time step. | |
213 | * \return dt. | |
214 | */ | |
215 | time_type current_time_step( void ) const | |
216 | { | |
217 | return m_dt; | |
218 | } | |
219 | ||
220 | ||
221 | private: | |
222 | ||
223 | state_type& get_current_state( void ) | |
224 | { | |
225 | return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ; | |
226 | } | |
227 | ||
228 | const state_type& get_current_state( void ) const | |
229 | { | |
230 | return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ; | |
231 | } | |
232 | ||
233 | state_type& get_old_state( void ) | |
234 | { | |
235 | return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ; | |
236 | } | |
237 | ||
238 | const state_type& get_old_state( void ) const | |
239 | { | |
240 | return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ; | |
241 | } | |
242 | ||
243 | void toggle_current_state( void ) | |
244 | { | |
245 | m_current_state_x1 = ! m_current_state_x1; | |
246 | } | |
247 | ||
248 | ||
249 | template< class StateIn > | |
250 | bool resize_impl( const StateIn &x ) | |
251 | { | |
252 | bool resized = false; | |
253 | resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() ); | |
254 | resized |= adjust_size_by_resizeability( m_x2 , x , typename is_resizeable<state_type>::type() ); | |
255 | return resized; | |
256 | } | |
257 | ||
258 | ||
259 | stepper_type m_stepper; | |
260 | resizer_type m_resizer; | |
261 | wrapped_state_type m_x1 , m_x2; | |
262 | bool m_current_state_x1; // if true, the current state is m_x1 | |
263 | time_type m_t , m_t_old , m_dt; | |
264 | ||
265 | }; | |
266 | ||
267 | ||
268 | ||
269 | ||
270 | ||
271 | /** | |
272 | * \brief The class representing dense-output Runge-Kutta steppers with FSAL property. | |
273 | * | |
274 | * The interface is the same as for dense_output_runge_kutta< Stepper , stepper_tag >. | |
275 | * This class provides dense output functionality based on methods with step size controlled | |
276 | * | |
277 | * | |
278 | * \tparam Stepper The stepper type of the underlying algorithm. | |
279 | */ | |
280 | template< class Stepper > | |
281 | class dense_output_runge_kutta< Stepper , explicit_controlled_stepper_fsal_tag > | |
282 | { | |
283 | public: | |
284 | ||
285 | /* | |
286 | * We do not need all typedefs. | |
287 | */ | |
288 | typedef Stepper controlled_stepper_type; | |
289 | ||
290 | typedef typename controlled_stepper_type::stepper_type stepper_type; | |
291 | typedef typename stepper_type::state_type state_type; | |
292 | typedef typename stepper_type::wrapped_state_type wrapped_state_type; | |
293 | typedef typename stepper_type::value_type value_type; | |
294 | typedef typename stepper_type::deriv_type deriv_type; | |
295 | typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type; | |
296 | typedef typename stepper_type::time_type time_type; | |
297 | typedef typename stepper_type::algebra_type algebra_type; | |
298 | typedef typename stepper_type::operations_type operations_type; | |
299 | typedef typename stepper_type::resizer_type resizer_type; | |
300 | typedef dense_output_stepper_tag stepper_category; | |
301 | typedef dense_output_runge_kutta< Stepper > dense_output_stepper_type; | |
302 | ||
303 | ||
304 | dense_output_runge_kutta( const controlled_stepper_type &stepper = controlled_stepper_type() ) | |
305 | : m_stepper( stepper ) , m_resizer() , | |
306 | m_current_state_x1( true ) , | |
307 | m_x1() , m_x2() , m_dxdt1() , m_dxdt2() , | |
308 | m_t() , m_t_old() , m_dt() , | |
309 | m_is_deriv_initialized( false ) | |
310 | { } | |
311 | ||
312 | ||
313 | template< class StateType > | |
314 | void initialize( const StateType &x0 , time_type t0 , time_type dt0 ) | |
315 | { | |
316 | m_resizer.adjust_size( x0 , detail::bind( &dense_output_stepper_type::template resize< StateType > , detail::ref( *this ) , detail::_1 ) ); | |
317 | boost::numeric::odeint::copy( x0 , get_current_state() ); | |
318 | m_t = t0; | |
319 | m_dt = dt0; | |
320 | m_is_deriv_initialized = false; | |
321 | } | |
322 | ||
323 | template< class System > | |
324 | std::pair< time_type , time_type > do_step( System system ) | |
325 | { | |
326 | if( !m_is_deriv_initialized ) | |
327 | { | |
328 | typename odeint::unwrap_reference< System >::type &sys = system; | |
329 | sys( get_current_state() , get_current_deriv() , m_t ); | |
330 | m_is_deriv_initialized = true; | |
331 | } | |
332 | ||
333 | failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails | |
334 | controlled_step_result res = fail; | |
335 | m_t_old = m_t; | |
336 | do | |
337 | { | |
338 | res = m_stepper.try_step( system , get_current_state() , get_current_deriv() , m_t , | |
339 | get_old_state() , get_old_deriv() , m_dt ); | |
340 | fail_checker(); // check for overflow of failed steps | |
341 | } | |
342 | while( res == fail ); | |
343 | toggle_current_state(); | |
344 | return std::make_pair( m_t_old , m_t ); | |
345 | } | |
346 | ||
347 | ||
348 | /* | |
349 | * The two overloads are needed in order to solve the forwarding problem. | |
350 | */ | |
351 | template< class StateOut > | |
352 | void calc_state( time_type t , StateOut &x ) const | |
353 | { | |
354 | m_stepper.stepper().calc_state( t , x , get_old_state() , get_old_deriv() , m_t_old , | |
355 | get_current_state() , get_current_deriv() , m_t ); | |
356 | } | |
357 | ||
358 | template< class StateOut > | |
359 | void calc_state( time_type t , const StateOut &x ) const | |
360 | { | |
361 | m_stepper.stepper().calc_state( t , x , get_old_state() , get_old_deriv() , m_t_old , | |
362 | get_current_state() , get_current_deriv() , m_t ); | |
363 | } | |
364 | ||
365 | ||
366 | template< class StateIn > | |
367 | bool resize( const StateIn &x ) | |
368 | { | |
369 | bool resized = false; | |
370 | resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() ); | |
371 | resized |= adjust_size_by_resizeability( m_x2 , x , typename is_resizeable<state_type>::type() ); | |
372 | resized |= adjust_size_by_resizeability( m_dxdt1 , x , typename is_resizeable<deriv_type>::type() ); | |
373 | resized |= adjust_size_by_resizeability( m_dxdt2 , x , typename is_resizeable<deriv_type>::type() ); | |
374 | return resized; | |
375 | } | |
376 | ||
377 | ||
378 | template< class StateType > | |
379 | void adjust_size( const StateType &x ) | |
380 | { | |
381 | resize( x ); | |
382 | m_stepper.stepper().resize( x ); | |
383 | } | |
384 | ||
385 | const state_type& current_state( void ) const | |
386 | { | |
387 | return get_current_state(); | |
388 | } | |
389 | ||
390 | time_type current_time( void ) const | |
391 | { | |
392 | return m_t; | |
393 | } | |
394 | ||
395 | const state_type& previous_state( void ) const | |
396 | { | |
397 | return get_old_state(); | |
398 | } | |
399 | ||
400 | time_type previous_time( void ) const | |
401 | { | |
402 | return m_t_old; | |
403 | } | |
404 | ||
405 | time_type current_time_step( void ) const | |
406 | { | |
407 | return m_dt; | |
408 | } | |
409 | ||
410 | ||
411 | private: | |
412 | ||
413 | state_type& get_current_state( void ) | |
414 | { | |
415 | return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ; | |
416 | } | |
417 | ||
418 | const state_type& get_current_state( void ) const | |
419 | { | |
420 | return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ; | |
421 | } | |
422 | ||
423 | state_type& get_old_state( void ) | |
424 | { | |
425 | return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ; | |
426 | } | |
427 | ||
428 | const state_type& get_old_state( void ) const | |
429 | { | |
430 | return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ; | |
431 | } | |
432 | ||
433 | deriv_type& get_current_deriv( void ) | |
434 | { | |
435 | return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ; | |
436 | } | |
437 | ||
438 | const deriv_type& get_current_deriv( void ) const | |
439 | { | |
440 | return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ; | |
441 | } | |
442 | ||
443 | deriv_type& get_old_deriv( void ) | |
444 | { | |
445 | return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ; | |
446 | } | |
447 | ||
448 | const deriv_type& get_old_deriv( void ) const | |
449 | { | |
450 | return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ; | |
451 | } | |
452 | ||
453 | ||
454 | void toggle_current_state( void ) | |
455 | { | |
456 | m_current_state_x1 = ! m_current_state_x1; | |
457 | } | |
458 | ||
459 | ||
460 | controlled_stepper_type m_stepper; | |
461 | resizer_type m_resizer; | |
462 | bool m_current_state_x1; | |
463 | wrapped_state_type m_x1 , m_x2; | |
464 | wrapped_deriv_type m_dxdt1 , m_dxdt2; | |
465 | time_type m_t , m_t_old , m_dt; | |
466 | bool m_is_deriv_initialized; | |
467 | ||
468 | }; | |
469 | ||
470 | } // namespace odeint | |
471 | } // namespace numeric | |
472 | } // namespace boost | |
473 | ||
474 | ||
475 | ||
476 | #endif // BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED |