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