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 */ dense_output_runge_kutta(const stepper_type & stepper=stepper_type ())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 > initialize(const StateType & x0,time_type t0,time_type dt0)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 > do_step(System 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(); 129 return std::make_pair( m_t_old , m_dt ); 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 > calc_state(time_type t,StateOut & x) const143 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 > calc_state(time_type t,const StateOut & x) const159 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 > adjust_size(const StateType & x)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 */ current_state(void) const179 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 */ current_time(void) const188 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 */ previous_state(void) const197 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 */ previous_time(void) const206 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 */ current_time_step(void) const215 time_type current_time_step( void ) const 216 { 217 return m_dt; 218 } 219 220 221 private: 222 get_current_state(void)223 state_type& get_current_state( void ) 224 { 225 return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ; 226 } 227 get_current_state(void) const228 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 get_old_state(void)233 state_type& get_old_state( void ) 234 { 235 return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ; 236 } 237 get_old_state(void) const238 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 toggle_current_state(void)243 void toggle_current_state( void ) 244 { 245 m_current_state_x1 = ! m_current_state_x1; 246 } 247 248 249 template< class StateIn > resize_impl(const StateIn & x)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 dense_output_runge_kutta(const controlled_stepper_type & stepper=controlled_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 ) 310 { } 311 312 313 template< class StateType > initialize(const StateType & x0,time_type t0,time_type dt0)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 > do_step(System 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 > calc_state(time_type t,StateOut & x) const352 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 > calc_state(time_type t,const StateOut & x) const359 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 > resize(const StateIn & x)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 > adjust_size(const StateType & x)379 void adjust_size( const StateType &x ) 380 { 381 resize( x ); 382 m_stepper.stepper().resize( x ); 383 } 384 current_state(void) const385 const state_type& current_state( void ) const 386 { 387 return get_current_state(); 388 } 389 current_time(void) const390 time_type current_time( void ) const 391 { 392 return m_t; 393 } 394 previous_state(void) const395 const state_type& previous_state( void ) const 396 { 397 return get_old_state(); 398 } 399 previous_time(void) const400 time_type previous_time( void ) const 401 { 402 return m_t_old; 403 } 404 current_time_step(void) const405 time_type current_time_step( void ) const 406 { 407 return m_dt; 408 } 409 410 411 private: 412 get_current_state(void)413 state_type& get_current_state( void ) 414 { 415 return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ; 416 } 417 get_current_state(void) const418 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 get_old_state(void)423 state_type& get_old_state( void ) 424 { 425 return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ; 426 } 427 get_old_state(void) const428 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 get_current_deriv(void)433 deriv_type& get_current_deriv( void ) 434 { 435 return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ; 436 } 437 get_current_deriv(void) const438 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 get_old_deriv(void)443 deriv_type& get_old_deriv( void ) 444 { 445 return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ; 446 } 447 get_old_deriv(void) const448 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 toggle_current_state(void)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 477