Chris@16
|
1 /*
|
Chris@16
|
2 [auto_generated]
|
Chris@16
|
3 boost/numeric/odeint/integrate/detail/integrate_const.hpp
|
Chris@16
|
4
|
Chris@16
|
5 [begin_description]
|
Chris@16
|
6 integrate const implementation
|
Chris@16
|
7 [end_description]
|
Chris@16
|
8
|
Chris@101
|
9 Copyright 2012 Mario Mulansky
|
Chris@101
|
10 Copyright 2012 Christoph Koke
|
Chris@101
|
11 Copyright 2012 Karsten Ahnert
|
Chris@16
|
12
|
Chris@16
|
13 Distributed under the Boost Software License, Version 1.0.
|
Chris@16
|
14 (See accompanying file LICENSE_1_0.txt or
|
Chris@16
|
15 copy at http://www.boost.org/LICENSE_1_0.txt)
|
Chris@16
|
16 */
|
Chris@16
|
17
|
Chris@16
|
18 #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
|
Chris@16
|
19 #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
|
Chris@16
|
20
|
Chris@16
|
21 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
|
Chris@16
|
22 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
|
Chris@16
|
23 #include <boost/numeric/odeint/util/unit_helper.hpp>
|
Chris@16
|
24 #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
|
Chris@16
|
25
|
Chris@16
|
26 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
|
Chris@16
|
27
|
Chris@16
|
28 namespace boost {
|
Chris@16
|
29 namespace numeric {
|
Chris@16
|
30 namespace odeint {
|
Chris@16
|
31 namespace detail {
|
Chris@16
|
32
|
Chris@16
|
33 // forward declaration
|
Chris@16
|
34 template< class Stepper , class System , class State , class Time , class Observer >
|
Chris@16
|
35 size_t integrate_adaptive(
|
Chris@16
|
36 Stepper stepper , System system , State &start_state ,
|
Chris@16
|
37 Time &start_time , Time end_time , Time &dt ,
|
Chris@16
|
38 Observer observer , controlled_stepper_tag
|
Chris@16
|
39 );
|
Chris@16
|
40
|
Chris@16
|
41
|
Chris@16
|
42 template< class Stepper , class System , class State , class Time , class Observer >
|
Chris@16
|
43 size_t integrate_const(
|
Chris@16
|
44 Stepper stepper , System system , State &start_state ,
|
Chris@16
|
45 Time start_time , Time end_time , Time dt ,
|
Chris@16
|
46 Observer observer , stepper_tag
|
Chris@16
|
47 )
|
Chris@16
|
48 {
|
Chris@16
|
49 typename odeint::unwrap_reference< Observer >::type &obs = observer;
|
Chris@101
|
50 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
|
Chris@16
|
51
|
Chris@16
|
52 Time time = start_time;
|
Chris@16
|
53 int step = 0;
|
Chris@101
|
54 // cast time+dt explicitely in case of expression templates (e.g. multiprecision)
|
Chris@101
|
55 while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
|
Chris@16
|
56 {
|
Chris@16
|
57 obs( start_state , time );
|
Chris@101
|
58 st.do_step( system , start_state , time , dt );
|
Chris@16
|
59 // direct computation of the time avoids error propagation happening when using time += dt
|
Chris@16
|
60 // we need clumsy type analysis to get boost units working here
|
Chris@16
|
61 ++step;
|
Chris@16
|
62 time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * dt;
|
Chris@16
|
63 }
|
Chris@16
|
64 obs( start_state , time );
|
Chris@16
|
65
|
Chris@16
|
66 return step;
|
Chris@16
|
67 }
|
Chris@16
|
68
|
Chris@16
|
69
|
Chris@16
|
70
|
Chris@16
|
71 template< class Stepper , class System , class State , class Time , class Observer >
|
Chris@16
|
72 size_t integrate_const(
|
Chris@16
|
73 Stepper stepper , System system , State &start_state ,
|
Chris@16
|
74 Time start_time , Time end_time , Time dt ,
|
Chris@16
|
75 Observer observer , controlled_stepper_tag
|
Chris@16
|
76 )
|
Chris@16
|
77 {
|
Chris@16
|
78 typename odeint::unwrap_reference< Observer >::type &obs = observer;
|
Chris@16
|
79
|
Chris@16
|
80 Time time = start_time;
|
Chris@16
|
81 const Time time_step = dt;
|
Chris@101
|
82 int real_steps = 0;
|
Chris@16
|
83 int step = 0;
|
Chris@16
|
84
|
Chris@101
|
85 while( less_eq_with_sign( static_cast<Time>(time+time_step) , end_time , dt ) )
|
Chris@16
|
86 {
|
Chris@16
|
87 obs( start_state , time );
|
Chris@101
|
88 real_steps += detail::integrate_adaptive( stepper , system , start_state , time , time+time_step , dt ,
|
Chris@16
|
89 null_observer() , controlled_stepper_tag() );
|
Chris@16
|
90 // direct computation of the time avoids error propagation happening when using time += dt
|
Chris@16
|
91 // we need clumsy type analysis to get boost units working here
|
Chris@101
|
92 step++;
|
Chris@16
|
93 time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * time_step;
|
Chris@16
|
94 }
|
Chris@16
|
95 obs( start_state , time );
|
Chris@101
|
96
|
Chris@101
|
97 return real_steps;
|
Chris@16
|
98 }
|
Chris@16
|
99
|
Chris@16
|
100
|
Chris@16
|
101 template< class Stepper , class System , class State , class Time , class Observer >
|
Chris@16
|
102 size_t integrate_const(
|
Chris@16
|
103 Stepper stepper , System system , State &start_state ,
|
Chris@16
|
104 Time start_time , Time end_time , Time dt ,
|
Chris@16
|
105 Observer observer , dense_output_stepper_tag
|
Chris@16
|
106 )
|
Chris@16
|
107 {
|
Chris@16
|
108 typename odeint::unwrap_reference< Observer >::type &obs = observer;
|
Chris@101
|
109 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
|
Chris@16
|
110
|
Chris@16
|
111 Time time = start_time;
|
Chris@16
|
112
|
Chris@101
|
113 st.initialize( start_state , time , dt );
|
Chris@16
|
114 obs( start_state , time );
|
Chris@16
|
115 time += dt;
|
Chris@16
|
116
|
Chris@16
|
117 int obs_step( 1 );
|
Chris@16
|
118 int real_step( 0 );
|
Chris@16
|
119
|
Chris@101
|
120 while( less_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
|
Chris@16
|
121 {
|
Chris@101
|
122 while( less_eq_with_sign( time , st.current_time() , dt ) )
|
Chris@16
|
123 {
|
Chris@101
|
124 st.calc_state( time , start_state );
|
Chris@16
|
125 obs( start_state , time );
|
Chris@16
|
126 ++obs_step;
|
Chris@16
|
127 // direct computation of the time avoids error propagation happening when using time += dt
|
Chris@16
|
128 // we need clumsy type analysis to get boost units working here
|
Chris@16
|
129 time = start_time + static_cast< typename unit_value_type<Time>::type >(obs_step) * dt;
|
Chris@16
|
130 }
|
Chris@16
|
131 // we have not reached the end, do another real step
|
Chris@101
|
132 if( less_with_sign( static_cast<Time>(st.current_time()+st.current_time_step()) ,
|
Chris@16
|
133 end_time ,
|
Chris@101
|
134 st.current_time_step() ) )
|
Chris@16
|
135 {
|
Chris@101
|
136 while( less_eq_with_sign( st.current_time() , time , dt ) )
|
Chris@16
|
137 {
|
Chris@101
|
138 st.do_step( system );
|
Chris@16
|
139 ++real_step;
|
Chris@16
|
140 }
|
Chris@16
|
141 }
|
Chris@101
|
142 else if( less_with_sign( st.current_time() , end_time , st.current_time_step() ) )
|
Chris@16
|
143 { // do the last step ending exactly on the end point
|
Chris@101
|
144 st.initialize( st.current_state() , st.current_time() , end_time - st.current_time() );
|
Chris@101
|
145 st.do_step( system );
|
Chris@16
|
146 ++real_step;
|
Chris@16
|
147 }
|
Chris@16
|
148
|
Chris@16
|
149 }
|
Chris@16
|
150 // last observation, if we are still in observation interval
|
Chris@16
|
151 if( less_eq_with_sign( time , end_time , dt ) )
|
Chris@16
|
152 {
|
Chris@101
|
153 st.calc_state( time , start_state );
|
Chris@16
|
154 obs( start_state , time );
|
Chris@16
|
155 }
|
Chris@16
|
156
|
Chris@16
|
157 return real_step;
|
Chris@16
|
158 }
|
Chris@16
|
159
|
Chris@16
|
160
|
Chris@16
|
161 } } } }
|
Chris@16
|
162
|
Chris@16
|
163 #endif
|