integrate_const.hpp 5.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159
  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 2009-2012 Karsten Ahnert
  8. Copyright 2009-2012 Mario Mulansky
  9. Distributed under the Boost Software License, Version 1.0.
  10. (See accompanying file LICENSE_1_0.txt or
  11. copy at http://www.boost.org/LICENSE_1_0.txt)
  12. */
  13. #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
  14. #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
  15. #include <boost/numeric/odeint/util/unwrap_reference.hpp>
  16. #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
  17. #include <boost/numeric/odeint/util/unit_helper.hpp>
  18. #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
  19. #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
  20. namespace boost {
  21. namespace numeric {
  22. namespace odeint {
  23. namespace detail {
  24. // forward declaration
  25. template< class Stepper , class System , class State , class Time , class Observer >
  26. size_t integrate_adaptive(
  27. Stepper stepper , System system , State &start_state ,
  28. Time &start_time , Time end_time , Time &dt ,
  29. Observer observer , controlled_stepper_tag
  30. );
  31. template< class Stepper , class System , class State , class Time , class Observer >
  32. size_t integrate_const(
  33. Stepper stepper , System system , State &start_state ,
  34. Time start_time , Time end_time , Time dt ,
  35. Observer observer , stepper_tag
  36. )
  37. {
  38. typename odeint::unwrap_reference< Observer >::type &obs = observer;
  39. Time time = start_time;
  40. int step = 0;
  41. while( less_eq_with_sign( time+dt , end_time , dt ) )
  42. {
  43. obs( start_state , time );
  44. stepper.do_step( system , start_state , time , dt );
  45. // direct computation of the time avoids error propagation happening when using time += dt
  46. // we need clumsy type analysis to get boost units working here
  47. ++step;
  48. time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * dt;
  49. }
  50. obs( start_state , time );
  51. return step;
  52. }
  53. template< class Stepper , class System , class State , class Time , class Observer >
  54. size_t integrate_const(
  55. Stepper stepper , System system , State &start_state ,
  56. Time start_time , Time end_time , Time dt ,
  57. Observer observer , controlled_stepper_tag
  58. )
  59. {
  60. typename odeint::unwrap_reference< Observer >::type &obs = observer;
  61. Time time = start_time;
  62. const Time time_step = dt;
  63. int step = 0;
  64. while( less_eq_with_sign( time+time_step , end_time , dt ) )
  65. {
  66. obs( start_state , time );
  67. detail::integrate_adaptive( stepper , system , start_state , 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. ++step;
  72. time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * time_step;
  73. }
  74. obs( start_state , time );
  75. return step;
  76. }
  77. template< class Stepper , class System , class State , class Time , class Observer >
  78. size_t integrate_const(
  79. Stepper stepper , System system , State &start_state ,
  80. Time start_time , Time end_time , Time dt ,
  81. Observer observer , dense_output_stepper_tag
  82. )
  83. {
  84. typename odeint::unwrap_reference< Observer >::type &obs = observer;
  85. Time time = start_time;
  86. stepper.initialize( start_state , time , dt );
  87. obs( start_state , time );
  88. time += dt;
  89. int obs_step( 1 );
  90. int real_step( 0 );
  91. while( less_with_sign( time+dt , end_time , dt ) )
  92. {
  93. while( less_eq_with_sign( time , stepper.current_time() , dt ) )
  94. {
  95. stepper.calc_state( time , start_state );
  96. obs( start_state , time );
  97. ++obs_step;
  98. // direct computation of the time avoids error propagation happening when using time += dt
  99. // we need clumsy type analysis to get boost units working here
  100. time = start_time + static_cast< typename unit_value_type<Time>::type >(obs_step) * dt;
  101. }
  102. // we have not reached the end, do another real step
  103. if( less_with_sign( stepper.current_time()+stepper.current_time_step() ,
  104. end_time ,
  105. stepper.current_time_step() ) )
  106. {
  107. while( less_eq_with_sign( stepper.current_time() , time , dt ) )
  108. {
  109. stepper.do_step( system );
  110. ++real_step;
  111. }
  112. }
  113. else if( less_with_sign( stepper.current_time() , end_time , stepper.current_time_step() ) )
  114. { // do the last step ending exactly on the end point
  115. stepper.initialize( stepper.current_state() , stepper.current_time() , end_time - stepper.current_time() );
  116. stepper.do_step( system );
  117. ++real_step;
  118. }
  119. }
  120. // last observation, if we are still in observation interval
  121. if( less_eq_with_sign( time , end_time , dt ) )
  122. {
  123. stepper.calc_state( time , start_state );
  124. obs( start_state , time );
  125. }
  126. return real_step;
  127. }
  128. } } } }
  129. #endif