Chris@16
|
1 /*
|
Chris@16
|
2 [auto_generated]
|
Chris@16
|
3 boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp
|
Chris@16
|
4
|
Chris@16
|
5 [begin_description]
|
Chris@16
|
6 Default Integrate adaptive implementation.
|
Chris@16
|
7 [end_description]
|
Chris@16
|
8
|
Chris@101
|
9 Copyright 2011-2013 Karsten Ahnert
|
Chris@101
|
10 Copyright 2011-2012 Mario Mulansky
|
Chris@101
|
11 Copyright 2012 Christoph Koke
|
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
|
Chris@16
|
19 #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED
|
Chris@16
|
20 #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED
|
Chris@16
|
21
|
Chris@16
|
22 #include <stdexcept>
|
Chris@16
|
23
|
Chris@101
|
24 #include <boost/throw_exception.hpp>
|
Chris@101
|
25
|
Chris@16
|
26 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
|
Chris@16
|
27 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
|
Chris@101
|
28 #include <boost/numeric/odeint/integrate/detail/integrate_const.hpp>
|
Chris@16
|
29 #include <boost/numeric/odeint/util/bind.hpp>
|
Chris@16
|
30 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
|
Chris@16
|
31 #include <boost/numeric/odeint/util/copy.hpp>
|
Chris@16
|
32
|
Chris@16
|
33 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
|
Chris@16
|
34
|
Chris@16
|
35
|
Chris@16
|
36 #include <iostream>
|
Chris@16
|
37
|
Chris@16
|
38 namespace boost {
|
Chris@16
|
39 namespace numeric {
|
Chris@16
|
40 namespace odeint {
|
Chris@16
|
41 namespace detail {
|
Chris@16
|
42
|
Chris@16
|
43 // forward declaration
|
Chris@16
|
44 template< class Stepper , class System , class State , class Time , class Observer>
|
Chris@101
|
45 size_t integrate_const(
|
Chris@16
|
46 Stepper stepper , System system , State &start_state ,
|
Chris@101
|
47 Time start_time , Time end_time , Time dt ,
|
Chris@16
|
48 Observer observer , stepper_tag );
|
Chris@16
|
49
|
Chris@16
|
50 /*
|
Chris@16
|
51 * integrate_adaptive for simple stepper is basically an integrate_const + some last step
|
Chris@16
|
52 */
|
Chris@16
|
53 template< class Stepper , class System , class State , class Time , class Observer >
|
Chris@16
|
54 size_t integrate_adaptive(
|
Chris@16
|
55 Stepper stepper , System system , State &start_state ,
|
Chris@16
|
56 Time start_time , Time end_time , Time dt ,
|
Chris@16
|
57 Observer observer , stepper_tag
|
Chris@16
|
58 )
|
Chris@16
|
59 {
|
Chris@101
|
60 size_t steps = detail::integrate_const( stepper , system , start_state , start_time ,
|
Chris@101
|
61 end_time , dt , observer , stepper_tag() );
|
Chris@101
|
62 typename odeint::unwrap_reference< Observer >::type &obs = observer;
|
Chris@101
|
63 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
|
Chris@101
|
64
|
Chris@101
|
65 Time end = start_time + dt*steps;
|
Chris@16
|
66 if( less_with_sign( end , end_time , dt ) )
|
Chris@16
|
67 { //make a last step to end exactly at end_time
|
Chris@101
|
68 st.do_step( system , start_state , end , end_time - end );
|
Chris@16
|
69 steps++;
|
Chris@16
|
70 obs( start_state , end_time );
|
Chris@16
|
71 }
|
Chris@16
|
72 return steps;
|
Chris@16
|
73 }
|
Chris@16
|
74
|
Chris@16
|
75
|
Chris@16
|
76 /*
|
Chris@16
|
77 * classical integrate adaptive
|
Chris@16
|
78 */
|
Chris@16
|
79 template< class Stepper , class System , class State , class Time , class Observer >
|
Chris@16
|
80 size_t integrate_adaptive(
|
Chris@16
|
81 Stepper stepper , System system , State &start_state ,
|
Chris@16
|
82 Time &start_time , Time end_time , Time &dt ,
|
Chris@16
|
83 Observer observer , controlled_stepper_tag
|
Chris@16
|
84 )
|
Chris@16
|
85 {
|
Chris@16
|
86 typename odeint::unwrap_reference< Observer >::type &obs = observer;
|
Chris@101
|
87 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
|
Chris@16
|
88
|
Chris@16
|
89 const size_t max_attempts = 1000;
|
Chris@16
|
90 const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found.";
|
Chris@16
|
91 size_t count = 0;
|
Chris@16
|
92 while( less_with_sign( start_time , end_time , dt ) )
|
Chris@16
|
93 {
|
Chris@16
|
94 obs( start_state , start_time );
|
Chris@101
|
95 if( less_with_sign( end_time , static_cast<Time>(start_time + dt) , dt ) )
|
Chris@16
|
96 {
|
Chris@16
|
97 dt = end_time - start_time;
|
Chris@16
|
98 }
|
Chris@16
|
99
|
Chris@16
|
100 size_t trials = 0;
|
Chris@16
|
101 controlled_step_result res = success;
|
Chris@16
|
102 do
|
Chris@16
|
103 {
|
Chris@101
|
104 res = st.try_step( system , start_state , start_time , dt );
|
Chris@16
|
105 ++trials;
|
Chris@16
|
106 }
|
Chris@16
|
107 while( ( res == fail ) && ( trials < max_attempts ) );
|
Chris@101
|
108 if( trials == max_attempts ) BOOST_THROW_EXCEPTION( std::overflow_error( error_string ) );
|
Chris@16
|
109
|
Chris@16
|
110 ++count;
|
Chris@16
|
111 }
|
Chris@16
|
112 obs( start_state , start_time );
|
Chris@16
|
113 return count;
|
Chris@16
|
114 }
|
Chris@16
|
115
|
Chris@16
|
116
|
Chris@16
|
117 /*
|
Chris@16
|
118 * integrate adaptive for dense output steppers
|
Chris@16
|
119 *
|
Chris@16
|
120 * step size control is used if the stepper supports it
|
Chris@16
|
121 */
|
Chris@16
|
122 template< class Stepper , class System , class State , class Time , class Observer >
|
Chris@16
|
123 size_t integrate_adaptive(
|
Chris@16
|
124 Stepper stepper , System system , State &start_state ,
|
Chris@16
|
125 Time start_time , Time end_time , Time dt ,
|
Chris@16
|
126 Observer observer , dense_output_stepper_tag )
|
Chris@16
|
127 {
|
Chris@16
|
128 typename odeint::unwrap_reference< Observer >::type &obs = observer;
|
Chris@101
|
129 typename odeint::unwrap_reference< Stepper >::type &st = stepper;
|
Chris@16
|
130
|
Chris@16
|
131 size_t count = 0;
|
Chris@101
|
132 st.initialize( start_state , start_time , dt );
|
Chris@16
|
133
|
Chris@101
|
134 while( less_with_sign( st.current_time() , end_time , st.current_time_step() ) )
|
Chris@16
|
135 {
|
Chris@101
|
136 while( less_eq_with_sign( static_cast<Time>(st.current_time() + st.current_time_step()) ,
|
Chris@16
|
137 end_time ,
|
Chris@101
|
138 st.current_time_step() ) )
|
Chris@16
|
139 { //make sure we don't go beyond the end_time
|
Chris@101
|
140 obs( st.current_state() , st.current_time() );
|
Chris@101
|
141 st.do_step( system );
|
Chris@16
|
142 ++count;
|
Chris@16
|
143 }
|
Chris@101
|
144 // calculate time step to arrive exactly at end time
|
Chris@101
|
145 st.initialize( st.current_state() , st.current_time() , static_cast<Time>(end_time - st.current_time()) );
|
Chris@16
|
146 }
|
Chris@101
|
147 obs( st.current_state() , st.current_time() );
|
Chris@16
|
148 // overwrite start_state with the final point
|
Chris@101
|
149 boost::numeric::odeint::copy( st.current_state() , start_state );
|
Chris@16
|
150 return count;
|
Chris@16
|
151 }
|
Chris@16
|
152
|
Chris@16
|
153
|
Chris@16
|
154
|
Chris@16
|
155
|
Chris@16
|
156 } // namespace detail
|
Chris@16
|
157 } // namespace odeint
|
Chris@16
|
158 } // namespace numeric
|
Chris@16
|
159 } // namespace boost
|
Chris@16
|
160
|
Chris@16
|
161
|
Chris@16
|
162 #endif // BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED
|