Chris@16
|
1 /*
|
Chris@16
|
2 [auto_generated]
|
Chris@16
|
3 boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
|
Chris@16
|
4
|
Chris@16
|
5 [begin_description]
|
Chris@16
|
6 Controller for the Rosenbrock4 method.
|
Chris@16
|
7 [end_description]
|
Chris@16
|
8
|
Chris@16
|
9 Copyright 2009-2011 Karsten Ahnert
|
Chris@16
|
10 Copyright 2009-2011 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
|
Chris@16
|
18 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED
|
Chris@16
|
19 #define BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED
|
Chris@16
|
20
|
Chris@16
|
21 #include <boost/config.hpp>
|
Chris@16
|
22 #include <boost/numeric/odeint/util/bind.hpp>
|
Chris@16
|
23
|
Chris@16
|
24 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
|
Chris@16
|
25 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
|
Chris@16
|
26
|
Chris@16
|
27 #include <boost/numeric/odeint/util/copy.hpp>
|
Chris@16
|
28 #include <boost/numeric/odeint/util/is_resizeable.hpp>
|
Chris@16
|
29
|
Chris@16
|
30 #include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
|
Chris@16
|
31
|
Chris@16
|
32 namespace boost {
|
Chris@16
|
33 namespace numeric {
|
Chris@16
|
34 namespace odeint {
|
Chris@16
|
35
|
Chris@16
|
36 template< class Stepper >
|
Chris@16
|
37 class rosenbrock4_controller
|
Chris@16
|
38 {
|
Chris@16
|
39 private:
|
Chris@16
|
40
|
Chris@16
|
41
|
Chris@16
|
42 public:
|
Chris@16
|
43
|
Chris@16
|
44 typedef Stepper stepper_type;
|
Chris@16
|
45 typedef typename stepper_type::value_type value_type;
|
Chris@16
|
46 typedef typename stepper_type::state_type state_type;
|
Chris@16
|
47 typedef typename stepper_type::wrapped_state_type wrapped_state_type;
|
Chris@16
|
48 typedef typename stepper_type::time_type time_type;
|
Chris@16
|
49 typedef typename stepper_type::deriv_type deriv_type;
|
Chris@16
|
50 typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
|
Chris@16
|
51 typedef typename stepper_type::resizer_type resizer_type;
|
Chris@16
|
52 typedef controlled_stepper_tag stepper_category;
|
Chris@16
|
53
|
Chris@16
|
54 typedef rosenbrock4_controller< Stepper > controller_type;
|
Chris@16
|
55
|
Chris@16
|
56
|
Chris@16
|
57 rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 , const stepper_type &stepper = stepper_type() )
|
Chris@16
|
58 : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) ,
|
Chris@16
|
59 m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
|
Chris@16
|
60 m_last_rejected( false )
|
Chris@16
|
61 { }
|
Chris@16
|
62
|
Chris@16
|
63
|
Chris@16
|
64 value_type error( const state_type &x , const state_type &xold , const state_type &xerr )
|
Chris@16
|
65 {
|
Chris@16
|
66 BOOST_USING_STD_MAX();
|
Chris@16
|
67 using std::abs;
|
Chris@16
|
68
|
Chris@16
|
69 const size_t n = x.size();
|
Chris@16
|
70 value_type err = 0.0 , sk = 0.0;
|
Chris@16
|
71 for( size_t i=0 ; i<n ; ++i )
|
Chris@16
|
72 {
|
Chris@16
|
73 sk = m_atol + m_rtol * max BOOST_PREVENT_MACRO_SUBSTITUTION ( abs( xold[i] ) , abs( x[i] ) );
|
Chris@16
|
74 err += xerr[i] * xerr[i] / sk / sk;
|
Chris@16
|
75 }
|
Chris@16
|
76 return sqrt( err / value_type( n ) );
|
Chris@16
|
77 }
|
Chris@16
|
78
|
Chris@16
|
79 value_type last_error( void ) const
|
Chris@16
|
80 {
|
Chris@16
|
81 return m_err_old;
|
Chris@16
|
82 }
|
Chris@16
|
83
|
Chris@16
|
84
|
Chris@16
|
85
|
Chris@16
|
86
|
Chris@16
|
87 template< class System >
|
Chris@16
|
88 boost::numeric::odeint::controlled_step_result
|
Chris@16
|
89 try_step( System sys , state_type &x , time_type &t , time_type &dt )
|
Chris@16
|
90 {
|
Chris@16
|
91 m_xnew_resizer.adjust_size( x , detail::bind( &controller_type::template resize_m_xnew< state_type > , detail::ref( *this ) , detail::_1 ) );
|
Chris@16
|
92 boost::numeric::odeint::controlled_step_result res = try_step( sys , x , t , m_xnew.m_v , dt );
|
Chris@16
|
93 if( res == success )
|
Chris@16
|
94 {
|
Chris@16
|
95 boost::numeric::odeint::copy( m_xnew.m_v , x );
|
Chris@16
|
96 }
|
Chris@16
|
97 return res;
|
Chris@16
|
98 }
|
Chris@16
|
99
|
Chris@16
|
100
|
Chris@16
|
101 template< class System >
|
Chris@16
|
102 boost::numeric::odeint::controlled_step_result
|
Chris@16
|
103 try_step( System sys , const state_type &x , time_type &t , state_type &xout , time_type &dt )
|
Chris@16
|
104 {
|
Chris@16
|
105 BOOST_USING_STD_MIN();
|
Chris@16
|
106 BOOST_USING_STD_MAX();
|
Chris@16
|
107 using std::pow;
|
Chris@16
|
108
|
Chris@16
|
109 static const value_type safe = 0.9 , fac1 = 5.0 , fac2 = 1.0 / 6.0;
|
Chris@16
|
110
|
Chris@16
|
111 m_xerr_resizer.adjust_size( x , detail::bind( &controller_type::template resize_m_xerr< state_type > , detail::ref( *this ) , detail::_1 ) );
|
Chris@16
|
112
|
Chris@16
|
113 m_stepper.do_step( sys , x , t , xout , dt , m_xerr.m_v );
|
Chris@16
|
114 value_type err = error( xout , x , m_xerr.m_v );
|
Chris@16
|
115
|
Chris@16
|
116 value_type fac = max BOOST_PREVENT_MACRO_SUBSTITUTION ( fac2 , min BOOST_PREVENT_MACRO_SUBSTITUTION ( fac1 , pow( err , 0.25 ) / safe ) );
|
Chris@16
|
117 value_type dt_new = dt / fac;
|
Chris@16
|
118 if ( err <= 1.0 )
|
Chris@16
|
119 {
|
Chris@16
|
120 if( m_first_step )
|
Chris@16
|
121 {
|
Chris@16
|
122 m_first_step = false;
|
Chris@16
|
123 }
|
Chris@16
|
124 else
|
Chris@16
|
125 {
|
Chris@16
|
126 value_type fac_pred = ( m_dt_old / dt ) * pow( err * err / m_err_old , 0.25 ) / safe;
|
Chris@16
|
127 fac_pred = max BOOST_PREVENT_MACRO_SUBSTITUTION ( fac2 , min BOOST_PREVENT_MACRO_SUBSTITUTION ( fac1 , fac_pred ) );
|
Chris@16
|
128 fac = max BOOST_PREVENT_MACRO_SUBSTITUTION ( fac , fac_pred );
|
Chris@16
|
129 dt_new = dt / fac;
|
Chris@16
|
130 }
|
Chris@16
|
131
|
Chris@16
|
132 m_dt_old = dt;
|
Chris@16
|
133 m_err_old = max BOOST_PREVENT_MACRO_SUBSTITUTION ( 0.01 , err );
|
Chris@16
|
134 if( m_last_rejected )
|
Chris@16
|
135 dt_new = ( dt >= 0.0 ? min BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) : max BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) );
|
Chris@16
|
136 t += dt;
|
Chris@16
|
137 dt = dt_new;
|
Chris@16
|
138 m_last_rejected = false;
|
Chris@16
|
139 return success;
|
Chris@16
|
140 }
|
Chris@16
|
141 else
|
Chris@16
|
142 {
|
Chris@16
|
143 dt = dt_new;
|
Chris@16
|
144 m_last_rejected = true;
|
Chris@16
|
145 return fail;
|
Chris@16
|
146 }
|
Chris@16
|
147 }
|
Chris@16
|
148
|
Chris@16
|
149
|
Chris@16
|
150 template< class StateType >
|
Chris@16
|
151 void adjust_size( const StateType &x )
|
Chris@16
|
152 {
|
Chris@16
|
153 resize_m_xerr( x );
|
Chris@16
|
154 resize_m_xnew( x );
|
Chris@16
|
155 }
|
Chris@16
|
156
|
Chris@16
|
157
|
Chris@16
|
158
|
Chris@16
|
159 stepper_type& stepper( void )
|
Chris@16
|
160 {
|
Chris@16
|
161 return m_stepper;
|
Chris@16
|
162 }
|
Chris@16
|
163
|
Chris@16
|
164 const stepper_type& stepper( void ) const
|
Chris@16
|
165 {
|
Chris@16
|
166 return m_stepper;
|
Chris@16
|
167 }
|
Chris@16
|
168
|
Chris@16
|
169
|
Chris@16
|
170
|
Chris@16
|
171
|
Chris@16
|
172 private:
|
Chris@16
|
173
|
Chris@16
|
174 template< class StateIn >
|
Chris@16
|
175 bool resize_m_xerr( const StateIn &x )
|
Chris@16
|
176 {
|
Chris@16
|
177 return adjust_size_by_resizeability( m_xerr , x , typename is_resizeable<state_type>::type() );
|
Chris@16
|
178 }
|
Chris@16
|
179
|
Chris@16
|
180 template< class StateIn >
|
Chris@16
|
181 bool resize_m_xnew( const StateIn &x )
|
Chris@16
|
182 {
|
Chris@16
|
183 return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
|
Chris@16
|
184 }
|
Chris@16
|
185
|
Chris@16
|
186
|
Chris@16
|
187 stepper_type m_stepper;
|
Chris@16
|
188 resizer_type m_xerr_resizer;
|
Chris@16
|
189 resizer_type m_xnew_resizer;
|
Chris@16
|
190 wrapped_state_type m_xerr;
|
Chris@16
|
191 wrapped_state_type m_xnew;
|
Chris@16
|
192 value_type m_atol , m_rtol;
|
Chris@16
|
193 bool m_first_step;
|
Chris@16
|
194 value_type m_err_old , m_dt_old;
|
Chris@16
|
195 bool m_last_rejected;
|
Chris@16
|
196 };
|
Chris@16
|
197
|
Chris@16
|
198
|
Chris@16
|
199
|
Chris@16
|
200
|
Chris@16
|
201
|
Chris@16
|
202
|
Chris@16
|
203 } // namespace odeint
|
Chris@16
|
204 } // namespace numeric
|
Chris@16
|
205 } // namespace boost
|
Chris@16
|
206
|
Chris@16
|
207
|
Chris@16
|
208 #endif // BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_INCLUDED
|