3 boost/numeric/odeint/integrate/detail/integrate_const.hpp
6 integrate const implementation
9 Copyright 2012-2015 Mario Mulansky
10 Copyright 2012 Christoph Koke
11 Copyright 2012 Karsten Ahnert
13 Distributed under the Boost Software License, Version 1.0.
14 (See accompanying file LICENSE_1_0.txt or
15 copy at http://www.boost.org/LICENSE_1_0.txt)
18 #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
19 #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
21 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
22 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
23 #include <boost/numeric/odeint/util/unit_helper.hpp>
24 #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
26 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
33 // forward declaration
34 template< class Stepper , class System , class State , class Time , class Observer >
35 size_t integrate_adaptive(
36 Stepper stepper , System system , State &start_state ,
37 Time &start_time , Time end_time , Time &dt ,
38 Observer observer , controlled_stepper_tag
42 template< class Stepper , class System , class State , class Time , class Observer >
43 size_t integrate_const(
44 Stepper stepper , System system , State &start_state ,
45 Time start_time , Time end_time , Time dt ,
46 Observer observer , stepper_tag
50 typename odeint::unwrap_reference< Observer >::type &obs = observer;
51 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
53 Time time = start_time;
55 // cast time+dt explicitely in case of expression templates (e.g. multiprecision)
56 while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
58 obs( start_state , time );
59 st.do_step( system , start_state , time , dt );
60 // direct computation of the time avoids error propagation happening when using time += dt
61 // we need clumsy type analysis to get boost units working here
63 time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * dt;
65 obs( start_state , time );
72 template< class Stepper , class System , class State , class Time , class Observer >
73 size_t integrate_const(
74 Stepper stepper , System system , State &start_state ,
75 Time start_time , Time end_time , Time dt ,
76 Observer observer , controlled_stepper_tag
79 typename odeint::unwrap_reference< Observer >::type &obs = observer;
81 Time time = start_time;
82 const Time time_step = dt;
86 while( less_eq_with_sign( static_cast<Time>(time+time_step) , end_time , dt ) )
88 obs( start_state , time );
89 // integrate_adaptive_checked uses the given checker to throw if an overflow occurs
90 real_steps += detail::integrate_adaptive(stepper, system, start_state, time,
91 static_cast<Time>(time + time_step), dt,
92 null_observer(), controlled_stepper_tag());
93 // direct computation of the time avoids error propagation happening when using time += dt
94 // we need clumsy type analysis to get boost units working here
96 time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * time_step;
98 obs( start_state , time );
104 template< class Stepper , class System , class State , class Time , class Observer >
105 size_t integrate_const(
106 Stepper stepper , System system , State &start_state ,
107 Time start_time , Time end_time , Time dt ,
108 Observer observer , dense_output_stepper_tag
111 typename odeint::unwrap_reference< Observer >::type &obs = observer;
112 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
114 Time time = start_time;
116 st.initialize( start_state , time , dt );
117 obs( start_state , time );
123 while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
125 while( less_eq_with_sign( time , st.current_time() , dt ) )
127 st.calc_state( time , start_state );
128 obs( start_state , time );
130 // direct computation of the time avoids error propagation happening when using time += dt
131 // we need clumsy type analysis to get boost units working here
132 time = start_time + static_cast< typename unit_value_type<Time>::type >(obs_step) * dt;
134 // we have not reached the end, do another real step
135 if( less_with_sign( static_cast<Time>(st.current_time()+st.current_time_step()) ,
137 st.current_time_step() ) )
139 while( less_eq_with_sign( st.current_time() , time , dt ) )
141 st.do_step( system );
145 else if( less_with_sign( st.current_time() , end_time , st.current_time_step() ) )
146 { // do the last step ending exactly on the end point
147 st.initialize( st.current_state() , st.current_time() , end_time - st.current_time() );
148 st.do_step( system );
153 // last observation, if we are still in observation interval
154 // might happen due to finite precision problems when computing the the time
155 if( less_eq_with_sign( time , end_time , dt ) )
157 st.calc_state( time , start_state );
158 obs( start_state , time );