integrate_n_steps.hpp 5.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161
  1. /*
  2. [auto_generated]
  3. boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp
  4. [begin_description]
  5. integrate steps implementation
  6. [end_description]
  7. Copyright 2012-2015 Mario Mulansky
  8. Copyright 2012 Christoph Koke
  9. Copyright 2012 Karsten Ahnert
  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_N_STEPS_HPP_INCLUDED
  15. #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_N_STEPS_HPP_INCLUDED
  16. #include <boost/numeric/odeint/util/unwrap_reference.hpp>
  17. #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
  18. #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
  19. #include <boost/numeric/odeint/util/unit_helper.hpp>
  20. #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
  21. namespace boost {
  22. namespace numeric {
  23. namespace odeint {
  24. namespace detail {
  25. // forward declaration
  26. template< class Stepper , class System , class State , class Time , class Observer >
  27. size_t integrate_adaptive_checked(
  28. Stepper stepper , System system , State &start_state ,
  29. Time &start_time , Time end_time , Time &dt ,
  30. Observer observer, controlled_stepper_tag
  31. );
  32. /* basic version */
  33. template< class Stepper , class System , class State , class Time , class Observer>
  34. Time integrate_n_steps(
  35. Stepper stepper , System system , State &start_state ,
  36. Time start_time , Time dt , size_t num_of_steps ,
  37. Observer observer , stepper_tag )
  38. {
  39. typename odeint::unwrap_reference< Observer >::type &obs = observer;
  40. typename odeint::unwrap_reference< Stepper >::type &st = stepper;
  41. Time time = start_time;
  42. for( size_t step = 0; step < num_of_steps ; ++step )
  43. {
  44. obs( start_state , time );
  45. st.do_step( system , start_state , time , dt );
  46. // direct computation of the time avoids error propagation happening when using time += dt
  47. // we need clumsy type analysis to get boost units working here
  48. time = start_time + static_cast< typename unit_value_type<Time>::type >( step+1 ) * dt;
  49. }
  50. obs( start_state , time );
  51. return time;
  52. }
  53. /* controlled version */
  54. template< class Stepper , class System , class State , class Time , class Observer >
  55. Time integrate_n_steps(
  56. Stepper stepper , System system , State &start_state ,
  57. Time start_time , Time dt , size_t num_of_steps ,
  58. Observer observer , controlled_stepper_tag )
  59. {
  60. typename odeint::unwrap_reference< Observer >::type &obs = observer;
  61. Time time = start_time;
  62. Time time_step = dt;
  63. for( size_t step = 0; step < num_of_steps ; ++step )
  64. {
  65. obs( start_state , time );
  66. // integrate_adaptive_checked uses the given checker to throw if an overflow occurs
  67. detail::integrate_adaptive(stepper, system, start_state, time, static_cast<Time>(time + time_step), dt,
  68. null_observer(), controlled_stepper_tag());
  69. // direct computation of the time avoids error propagation happening when using time += dt
  70. // we need clumsy type analysis to get boost units working here
  71. time = start_time + static_cast< typename unit_value_type<Time>::type >(step+1) * time_step;
  72. }
  73. obs( start_state , time );
  74. return time;
  75. }
  76. /* dense output version */
  77. template< class Stepper , class System , class State , class Time , class Observer >
  78. Time integrate_n_steps(
  79. Stepper stepper , System system , State &start_state ,
  80. Time start_time , Time dt , size_t num_of_steps ,
  81. Observer observer , dense_output_stepper_tag )
  82. {
  83. typename odeint::unwrap_reference< Observer >::type &obs = observer;
  84. typename odeint::unwrap_reference< Stepper >::type &st = stepper;
  85. Time time = start_time;
  86. const Time end_time = start_time + static_cast< typename unit_value_type<Time>::type >(num_of_steps) * dt;
  87. st.initialize( start_state , time , dt );
  88. size_t step = 0;
  89. while( step < num_of_steps )
  90. {
  91. while( less_with_sign( time , st.current_time() , st.current_time_step() ) )
  92. {
  93. st.calc_state( time , start_state );
  94. obs( start_state , time );
  95. ++step;
  96. // direct computation of the time avoids error propagation happening when using time += dt
  97. // we need clumsy type analysis to get boost units working here
  98. time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * dt;
  99. }
  100. // we have not reached the end, do another real step
  101. if( less_with_sign( static_cast<Time>(st.current_time()+st.current_time_step()) ,
  102. end_time ,
  103. st.current_time_step() ) )
  104. {
  105. st.do_step( system );
  106. }
  107. else if( less_with_sign( st.current_time() , end_time , st.current_time_step() ) )
  108. { // do the last step ending exactly on the end point
  109. st.initialize( st.current_state() , st.current_time() , static_cast<Time>(end_time - st.current_time()) );
  110. st.do_step( system );
  111. }
  112. }
  113. // make sure we really end exactly where we should end
  114. while( st.current_time() < end_time )
  115. {
  116. if( less_with_sign( end_time ,
  117. static_cast<Time>(st.current_time()+st.current_time_step()) ,
  118. st.current_time_step() ) )
  119. st.initialize( st.current_state() , st.current_time() , static_cast<Time>(end_time - st.current_time()) );
  120. st.do_step( system );
  121. }
  122. // observation at final point
  123. obs( st.current_state() , end_time );
  124. return time;
  125. }
  126. }
  127. }
  128. }
  129. }
  130. #endif /* BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_N_STEPS_HPP_INCLUDED */