integrate_const.hpp 5.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167
  1. /*
  2. [auto_generated]
  3. boost/numeric/odeint/integrate/detail/integrate_const.hpp
  4. [begin_description]
  5. integrate const 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_CONST_HPP_INCLUDED
  15. #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_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/util/unit_helper.hpp>
  19. #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.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(
  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. template< class Stepper , class System , class State , class Time , class Observer >
  33. size_t integrate_const(
  34. Stepper stepper , System system , State &start_state ,
  35. Time start_time , Time end_time , Time dt ,
  36. Observer observer , stepper_tag
  37. )
  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. int step = 0;
  43. // cast time+dt explicitely in case of expression templates (e.g. multiprecision)
  44. while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
  45. {
  46. obs( start_state , time );
  47. st.do_step( system , start_state , time , dt );
  48. // direct computation of the time avoids error propagation happening when using time += dt
  49. // we need clumsy type analysis to get boost units working here
  50. ++step;
  51. time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * dt;
  52. }
  53. obs( start_state , time );
  54. return step;
  55. }
  56. template< class Stepper , class System , class State , class Time , class Observer >
  57. size_t integrate_const(
  58. Stepper stepper , System system , State &start_state ,
  59. Time start_time , Time end_time , Time dt ,
  60. Observer observer , controlled_stepper_tag
  61. )
  62. {
  63. typename odeint::unwrap_reference< Observer >::type &obs = observer;
  64. Time time = start_time;
  65. const Time time_step = dt;
  66. int real_steps = 0;
  67. int step = 0;
  68. while( less_eq_with_sign( static_cast<Time>(time+time_step) , end_time , dt ) )
  69. {
  70. obs( start_state , time );
  71. // integrate_adaptive_checked uses the given checker to throw if an overflow occurs
  72. real_steps += detail::integrate_adaptive(stepper, system, start_state, time,
  73. static_cast<Time>(time + time_step), dt,
  74. null_observer(), controlled_stepper_tag());
  75. // direct computation of the time avoids error propagation happening when using time += dt
  76. // we need clumsy type analysis to get boost units working here
  77. step++;
  78. time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * time_step;
  79. }
  80. obs( start_state , time );
  81. return real_steps;
  82. }
  83. template< class Stepper , class System , class State , class Time , class Observer >
  84. size_t integrate_const(
  85. Stepper stepper , System system , State &start_state ,
  86. Time start_time , Time end_time , Time dt ,
  87. Observer observer , dense_output_stepper_tag
  88. )
  89. {
  90. typename odeint::unwrap_reference< Observer >::type &obs = observer;
  91. typename odeint::unwrap_reference< Stepper >::type &st = stepper;
  92. Time time = start_time;
  93. st.initialize( start_state , time , dt );
  94. obs( start_state , time );
  95. time += dt;
  96. int obs_step( 1 );
  97. int real_step( 0 );
  98. while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
  99. {
  100. while( less_eq_with_sign( time , st.current_time() , dt ) )
  101. {
  102. st.calc_state( time , start_state );
  103. obs( start_state , time );
  104. ++obs_step;
  105. // direct computation of the time avoids error propagation happening when using time += dt
  106. // we need clumsy type analysis to get boost units working here
  107. time = start_time + static_cast< typename unit_value_type<Time>::type >(obs_step) * dt;
  108. }
  109. // we have not reached the end, do another real step
  110. if( less_with_sign( static_cast<Time>(st.current_time()+st.current_time_step()) ,
  111. end_time ,
  112. st.current_time_step() ) )
  113. {
  114. while( less_eq_with_sign( st.current_time() , time , dt ) )
  115. {
  116. st.do_step( system );
  117. ++real_step;
  118. }
  119. }
  120. else if( less_with_sign( st.current_time() , end_time , st.current_time_step() ) )
  121. { // do the last step ending exactly on the end point
  122. st.initialize( st.current_state() , st.current_time() , end_time - st.current_time() );
  123. st.do_step( system );
  124. ++real_step;
  125. }
  126. }
  127. // last observation, if we are still in observation interval
  128. // might happen due to finite precision problems when computing the the time
  129. if( less_eq_with_sign( time , end_time , dt ) )
  130. {
  131. st.calc_state( time , start_state );
  132. obs( start_state , time );
  133. }
  134. return real_step;
  135. }
  136. } } } }
  137. #endif