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