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