integrate_times.hpp 5.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179
  1. /*
  2. [auto_generated]
  3. boost/numeric/odeint/integrate/detail/integrate_times.hpp
  4. [begin_description]
  5. Default integrate times implementation.
  6. [end_description]
  7. Copyright 2011-2015 Mario Mulansky
  8. Copyright 2012 Karsten Ahnert
  9. Copyright 2012 Christoph Koke
  10. Distributed under the Boost Software License, Version 1.0.
  11. (See accompanying file LICENSE_1_0.txt or
  12. copy at http://www.boost.org/LICENSE_1_0.txt)
  13. */
  14. #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
  15. #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
  16. #include <stdexcept>
  17. #include <boost/config.hpp>
  18. #include <boost/throw_exception.hpp>
  19. #include <boost/numeric/odeint/util/unwrap_reference.hpp>
  20. #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
  21. #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
  22. #include <boost/numeric/odeint/integrate/max_step_checker.hpp>
  23. namespace boost {
  24. namespace numeric {
  25. namespace odeint {
  26. namespace detail {
  27. /*
  28. * integrate_times for simple stepper
  29. */
  30. template<class Stepper, class System, class State, class TimeIterator, class Time, class Observer>
  31. size_t integrate_times(
  32. Stepper stepper , System system , State &start_state ,
  33. TimeIterator start_time , TimeIterator end_time , Time dt ,
  34. Observer observer , stepper_tag
  35. )
  36. {
  37. typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
  38. typedef typename odeint::unwrap_reference< Observer >::type observer_type;
  39. stepper_type &st = stepper;
  40. observer_type &obs = observer;
  41. typedef typename unit_value_type<Time>::type time_type;
  42. size_t steps = 0;
  43. Time current_dt = dt;
  44. while( true )
  45. {
  46. Time current_time = *start_time++;
  47. obs( start_state , current_time );
  48. if( start_time == end_time )
  49. break;
  50. while( less_with_sign( current_time , static_cast<time_type>(*start_time) , current_dt ) )
  51. {
  52. current_dt = min_abs( dt , *start_time - current_time );
  53. st.do_step( system , start_state , current_time , current_dt );
  54. current_time += current_dt;
  55. steps++;
  56. }
  57. }
  58. return steps;
  59. }
  60. /*
  61. * integrate_times for controlled stepper
  62. */
  63. template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
  64. size_t integrate_times(
  65. Stepper stepper , System system , State &start_state ,
  66. TimeIterator start_time , TimeIterator end_time , Time dt ,
  67. Observer observer , controlled_stepper_tag
  68. )
  69. {
  70. typename odeint::unwrap_reference< Observer >::type &obs = observer;
  71. typename odeint::unwrap_reference< Stepper >::type &st = stepper;
  72. typedef typename unit_value_type<Time>::type time_type;
  73. failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails
  74. size_t steps = 0;
  75. while( true )
  76. {
  77. Time current_time = *start_time++;
  78. obs( start_state , current_time );
  79. if( start_time == end_time )
  80. break;
  81. while( less_with_sign( current_time , static_cast<time_type>(*start_time) , dt ) )
  82. {
  83. // adjust stepsize to end up exactly at the observation point
  84. Time current_dt = min_abs( dt , *start_time - current_time );
  85. if( st.try_step( system , start_state , current_time , current_dt ) == success )
  86. {
  87. ++steps;
  88. // successful step -> reset the fail counter, see #173
  89. fail_checker.reset();
  90. // continue with the original step size if dt was reduced due to observation
  91. dt = max_abs( dt , current_dt );
  92. }
  93. else
  94. {
  95. fail_checker(); // check for possible overflow of failed steps in step size adjustment
  96. dt = current_dt;
  97. }
  98. }
  99. }
  100. return steps;
  101. }
  102. /*
  103. * integrate_times for dense output stepper
  104. */
  105. template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
  106. size_t integrate_times(
  107. Stepper stepper , System system , State &start_state ,
  108. TimeIterator start_time , TimeIterator end_time , Time dt ,
  109. Observer observer , dense_output_stepper_tag
  110. )
  111. {
  112. typename odeint::unwrap_reference< Observer >::type &obs = observer;
  113. typename odeint::unwrap_reference< Stepper >::type &st = stepper;
  114. typedef typename unit_value_type<Time>::type time_type;
  115. if( start_time == end_time )
  116. return 0;
  117. TimeIterator last_time_iterator = end_time;
  118. --last_time_iterator;
  119. Time last_time_point = static_cast<time_type>(*last_time_iterator);
  120. st.initialize( start_state , *start_time , dt );
  121. obs( start_state , *start_time++ );
  122. size_t count = 0;
  123. while( start_time != end_time )
  124. {
  125. while( ( start_time != end_time ) && less_eq_with_sign( static_cast<time_type>(*start_time) , st.current_time() , st.current_time_step() ) )
  126. {
  127. st.calc_state( *start_time , start_state );
  128. obs( start_state , *start_time );
  129. start_time++;
  130. }
  131. // we have not reached the end, do another real step
  132. if( less_eq_with_sign( st.current_time() + st.current_time_step() ,
  133. last_time_point ,
  134. st.current_time_step() ) )
  135. {
  136. st.do_step( system );
  137. ++count;
  138. }
  139. else if( start_time != end_time )
  140. { // do the last step ending exactly on the end point
  141. st.initialize( st.current_state() , st.current_time() , last_time_point - st.current_time() );
  142. st.do_step( system );
  143. ++count;
  144. }
  145. }
  146. return count;
  147. }
  148. } // namespace detail
  149. } // namespace odeint
  150. } // namespace numeric
  151. } // namespace boost
  152. #endif // BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED