3 boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp
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.
10 Copyright 2011-2013 Karsten Ahnert
11 Copyright 2011-2015 Mario Mulansky
12 Copyright 2012 Christoph Koke
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)
20 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED
21 #define BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED
27 #include <boost/throw_exception.hpp>
29 #include <boost/numeric/odeint/util/bind.hpp>
31 #include <boost/numeric/odeint/util/copy.hpp>
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>
37 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
38 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
40 #include <boost/numeric/odeint/integrate/max_step_checker.hpp>
46 template< class Stepper , class StepperCategory = typename Stepper::stepper_category >
47 class dense_output_runge_kutta;
51 * \brief The class representing dense-output Runge-Kutta steppers.
52 * \note In this stepper, the initialize method has to be called before using
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
61 * \tparam Stepper The stepper type of the underlying algorithm.
63 template< class Stepper >
64 class dense_output_runge_kutta< Stepper , stepper_tag >
70 * We do not need all typedefs.
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;
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.
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()
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.
105 template< class StateType >
106 void initialize( const StateType &x0 , time_type t0 , time_type dt0 )
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() );
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.
122 template< class System >
123 std::pair< time_type , time_type > do_step( System system )
125 m_stepper.do_step( system , get_current_state() , m_t , get_old_state() , m_dt );
128 toggle_current_state();
129 return std::make_pair( m_t_old , m_dt );
133 * The next two overloads are needed to solve the forwarding problem
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.
142 template< class StateOut >
143 void calc_state( time_type t , StateOut &x ) const
145 if( t == current_time() )
147 boost::numeric::odeint::copy( get_current_state() , x );
149 m_stepper.calc_state( x , t , get_old_state() , m_t_old , get_current_state() , m_t );
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.
158 template< class StateOut >
159 void calc_state( time_type t , const StateOut &x ) const
161 m_stepper.calc_state( x , t , get_old_state() , m_t_old , get_current_state() , m_t );
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.
168 template< class StateType >
169 void adjust_size( const StateType &x )
172 m_stepper.stepper().resize( x );
176 * \brief Returns the current state of the solution.
177 * \return The current state of the solution x(t).
179 const state_type& current_state( void ) const
181 return get_current_state();
185 * \brief Returns the current time of the solution.
186 * \return The current time of the solution t.
188 time_type current_time( void ) const
194 * \brief Returns the last state of the solution.
195 * \return The last state of the solution x(t-dt).
197 const state_type& previous_state( void ) const
199 return get_old_state();
203 * \brief Returns the last time of the solution.
204 * \return The last time of the solution t-dt.
206 time_type previous_time( void ) const
212 * \brief Returns the current time step.
215 time_type current_time_step( void ) const
223 state_type& get_current_state( void )
225 return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
228 const state_type& get_current_state( void ) const
230 return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
233 state_type& get_old_state( void )
235 return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
238 const state_type& get_old_state( void ) const
240 return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
243 void toggle_current_state( void )
245 m_current_state_x1 = ! m_current_state_x1;
249 template< class StateIn >
250 bool resize_impl( const StateIn &x )
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() );
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;
272 * \brief The class representing dense-output Runge-Kutta steppers with FSAL property.
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
278 * \tparam Stepper The stepper type of the underlying algorithm.
280 template< class Stepper >
281 class dense_output_runge_kutta< Stepper , explicit_controlled_stepper_fsal_tag >
286 * We do not need all typedefs.
288 typedef Stepper controlled_stepper_type;
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;
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 )
313 template< class StateType >
314 void initialize( const StateType &x0 , time_type t0 , time_type dt0 )
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() );
320 m_is_deriv_initialized = false;
323 template< class System >
324 std::pair< time_type , time_type > do_step( System system )
326 if( !m_is_deriv_initialized )
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;
333 failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails
334 controlled_step_result res = fail;
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
342 while( res == fail );
343 toggle_current_state();
344 return std::make_pair( m_t_old , m_t );
349 * The two overloads are needed in order to solve the forwarding problem.
351 template< class StateOut >
352 void calc_state( time_type t , StateOut &x ) const
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 );
358 template< class StateOut >
359 void calc_state( time_type t , const StateOut &x ) const
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 );
366 template< class StateIn >
367 bool resize( const StateIn &x )
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() );
378 template< class StateType >
379 void adjust_size( const StateType &x )
382 m_stepper.stepper().resize( x );
385 const state_type& current_state( void ) const
387 return get_current_state();
390 time_type current_time( void ) const
395 const state_type& previous_state( void ) const
397 return get_old_state();
400 time_type previous_time( void ) const
405 time_type current_time_step( void ) const
413 state_type& get_current_state( void )
415 return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
418 const state_type& get_current_state( void ) const
420 return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ;
423 state_type& get_old_state( void )
425 return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
428 const state_type& get_old_state( void ) const
430 return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ;
433 deriv_type& get_current_deriv( void )
435 return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
438 const deriv_type& get_current_deriv( void ) const
440 return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ;
443 deriv_type& get_old_deriv( void )
445 return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
448 const deriv_type& get_old_deriv( void ) const
450 return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ;
454 void toggle_current_state( void )
456 m_current_state_x1 = ! m_current_state_x1;
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;
470 } // namespace odeint
471 } // namespace numeric
476 #endif // BOOST_NUMERIC_ODEINT_STEPPER_DENSE_OUTPUT_RUNGE_KUTTA_HPP_INCLUDED