1 /* 2 [auto_generated] 3 boost/numeric/odeint/stepper/extrapolation_stepper.hpp 4 5 [begin_description] 6 extrapolation stepper 7 [end_description] 8 9 Copyright 2009-2015 Mario Mulansky 10 11 Distributed under the Boost Software License, Version 1.0. 12 (See accompanying file LICENSE_1_0.txt or 13 copy at http://www.boost.org/LICENSE_1_0.txt) 14 */ 15 16 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED 17 #define BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED 18 19 #include <iostream> 20 21 #include <algorithm> 22 23 #include <boost/config.hpp> // for min/max guidelines 24 #include <boost/static_assert.hpp> 25 26 #include <boost/numeric/odeint/util/bind.hpp> 27 #include <boost/numeric/odeint/util/unwrap_reference.hpp> 28 29 #include <boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp> 30 #include <boost/numeric/odeint/stepper/modified_midpoint.hpp> 31 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp> 32 #include <boost/numeric/odeint/algebra/range_algebra.hpp> 33 #include <boost/numeric/odeint/algebra/default_operations.hpp> 34 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> 35 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> 36 37 #include <boost/numeric/odeint/util/state_wrapper.hpp> 38 #include <boost/numeric/odeint/util/is_resizeable.hpp> 39 #include <boost/numeric/odeint/util/resizer.hpp> 40 #include <boost/numeric/odeint/util/unit_helper.hpp> 41 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp> 42 43 namespace boost 44 { 45 namespace numeric 46 { 47 namespace odeint 48 { 49 50 template < unsigned short Order, class State, class Value = double, 51 class Deriv = State, class Time = Value, 52 class Algebra = typename algebra_dispatcher< State >::algebra_type, 53 class Operations = 54 typename operations_dispatcher< State >::operations_type, 55 class Resizer = initially_resizer > 56 #ifndef DOXYGEN_SKIP 57 class extrapolation_stepper 58 : public explicit_error_stepper_base< 59 extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra, 60 Operations, Resizer >, 61 Order, Order, Order - 2, State, Value, Deriv, Time, Algebra, 62 Operations, Resizer > 63 #else 64 class extrapolation_stepper : public explicit_error_stepper_base 65 #endif 66 { 67 68 private: 69 // check for Order being odd 70 BOOST_STATIC_ASSERT_MSG( 71 ( ( Order % 2 ) == 0 ) && ( Order > 2 ), 72 "extrapolation_stepper requires even Order larger than 2" ); 73 74 public: 75 #ifndef DOXYGEN_SKIP 76 typedef explicit_error_stepper_base< 77 extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra, 78 Operations, Resizer >, 79 Order, Order, Order - 2, State, Value, Deriv, Time, Algebra, Operations, 80 Resizer > stepper_base_type; 81 #else 82 typedef explicit_error_stepper_base< extrapolation_stepper< ... >, ... > 83 stepper_base_type; 84 #endif 85 86 typedef typename stepper_base_type::state_type state_type; 87 typedef typename stepper_base_type::value_type value_type; 88 typedef typename stepper_base_type::deriv_type deriv_type; 89 typedef typename stepper_base_type::time_type time_type; 90 typedef typename stepper_base_type::algebra_type algebra_type; 91 typedef typename stepper_base_type::operations_type operations_type; 92 typedef typename stepper_base_type::resizer_type resizer_type; 93 94 #ifndef DOXYGEN_SKIP 95 typedef typename stepper_base_type::stepper_type stepper_type; 96 typedef typename stepper_base_type::wrapped_state_type wrapped_state_type; 97 typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type; 98 99 typedef std::vector< value_type > value_vector; 100 typedef std::vector< value_vector > value_matrix; 101 typedef std::vector< size_t > int_vector; 102 typedef std::vector< wrapped_state_type > state_table_type; 103 typedef modified_midpoint< state_type, value_type, deriv_type, time_type, 104 algebra_type, operations_type, 105 resizer_type > midpoint_stepper_type; 106 107 #endif // DOXYGEN_SKIP 108 109 typedef unsigned short order_type; 110 static const order_type order_value = stepper_base_type::order_value; 111 static const order_type stepper_order_value = 112 stepper_base_type::stepper_order_value; 113 static const order_type error_order_value = 114 stepper_base_type::error_order_value; 115 116 const static size_t m_k_max = ( order_value - 2 ) / 2; 117 extrapolation_stepper(const algebra_type & algebra=algebra_type ())118 extrapolation_stepper( const algebra_type &algebra = algebra_type() ) 119 : stepper_base_type( algebra ), m_interval_sequence( m_k_max + 1 ), 120 m_coeff( m_k_max + 1 ), m_table( m_k_max ) 121 { 122 for ( unsigned short i = 0; i < m_k_max + 1; i++ ) 123 { 124 m_interval_sequence[i] = 2 * ( i + 1 ); 125 m_coeff[i].resize( i ); 126 for ( size_t k = 0; k < i; ++k ) 127 { 128 const value_type r = 129 static_cast< value_type >( m_interval_sequence[i] ) / 130 static_cast< value_type >( m_interval_sequence[k] ); 131 m_coeff[i][k] = 132 static_cast< value_type >( 1 ) / 133 ( r * r - static_cast< value_type >( 134 1 ) ); // coefficients for extrapolation 135 } 136 } 137 } 138 139 template < class System, class StateIn, class DerivIn, class StateOut, 140 class Err > do_step_impl(System system,const StateIn & in,const DerivIn & dxdt,time_type t,StateOut & out,time_type dt,Err & xerr)141 void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt, 142 time_type t, StateOut &out, time_type dt, Err &xerr ) 143 { 144 // std::cout << "dt: " << dt << std::endl; 145 // normal step 146 do_step_impl( system, in, dxdt, t, out, dt ); 147 148 static const value_type val1( 1.0 ); 149 // additionally, perform the error calculation 150 stepper_base_type::m_algebra.for_each3( 151 xerr, out, m_table[0].m_v, 152 typename operations_type::template scale_sum2< 153 value_type, value_type >( val1, -val1 ) ); 154 } 155 156 template < class System, class StateInOut, class DerivIn, class Err > do_step_impl_io(System system,StateInOut & inout,const DerivIn & dxdt,time_type t,time_type dt,Err & xerr)157 void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt, 158 time_type t, time_type dt, Err &xerr ) 159 { 160 // normal step 161 do_step_impl_io( system, inout, dxdt, t, dt ); 162 163 static const value_type val1( 1.0 ); 164 // additionally, perform the error calculation 165 stepper_base_type::m_algebra.for_each3( 166 xerr, inout, m_table[0].m_v, 167 typename operations_type::template scale_sum2< 168 value_type, value_type >( val1, -val1 ) ); 169 } 170 171 template < class System, class StateIn, class DerivIn, class StateOut > do_step_impl(System system,const StateIn & in,const DerivIn & dxdt,time_type t,StateOut & out,time_type dt)172 void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt, 173 time_type t, StateOut &out, time_type dt ) 174 { 175 m_resizer.adjust_size( 176 in, detail::bind( &stepper_type::template resize_impl< StateIn >, 177 detail::ref( *this ), detail::_1 ) ); 178 size_t k = 0; 179 m_midpoint.set_steps( m_interval_sequence[k] ); 180 m_midpoint.do_step( system, in, dxdt, t, out, dt ); 181 for ( k = 1; k <= m_k_max; ++k ) 182 { 183 m_midpoint.set_steps( m_interval_sequence[k] ); 184 m_midpoint.do_step( system, in, dxdt, t, m_table[k - 1].m_v, dt ); 185 extrapolate( k, m_table, m_coeff, out ); 186 } 187 } 188 189 template < class System, class StateInOut, class DerivIn > do_step_impl_io(System system,StateInOut & inout,const DerivIn & dxdt,time_type t,time_type dt)190 void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt, 191 time_type t, time_type dt ) 192 { 193 // special care for inout 194 m_xout_resizer.adjust_size( 195 inout, 196 detail::bind( &stepper_type::template resize_m_xout< StateInOut >, 197 detail::ref( *this ), detail::_1 ) ); 198 do_step_impl( system, inout, dxdt, t, m_xout.m_v, dt ); 199 boost::numeric::odeint::copy( m_xout.m_v, inout ); 200 } 201 202 template < class System, class StateInOut, class DerivIn > do_step_dxdt_impl(System system,StateInOut & x,const DerivIn & dxdt,time_type t,time_type dt)203 void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt, 204 time_type t, time_type dt ) 205 { 206 do_step_impl_io( system , x , dxdt , t , dt ); 207 } 208 209 template < class System, class StateIn, class DerivIn, class StateOut > do_step_dxdt_impl(System system,const StateIn & in,const DerivIn & dxdt,time_type t,StateOut & out,time_type dt)210 void do_step_dxdt_impl( System system, const StateIn &in, 211 const DerivIn &dxdt, time_type t, StateOut &out, 212 time_type dt ) 213 { 214 do_step_impl( system , in , dxdt , t , out , dt ); 215 } 216 217 adjust_size(const StateIn & x)218 template < class StateIn > void adjust_size( const StateIn &x ) 219 { 220 resize_impl( x ); 221 m_midpoint.adjust_size( x ); 222 } 223 224 private: resize_impl(const StateIn & x)225 template < class StateIn > bool resize_impl( const StateIn &x ) 226 { 227 bool resized( false ); 228 for ( size_t i = 0; i < m_k_max; ++i ) 229 resized |= adjust_size_by_resizeability( 230 m_table[i], x, typename is_resizeable< state_type >::type() ); 231 return resized; 232 } 233 resize_m_xout(const StateIn & x)234 template < class StateIn > bool resize_m_xout( const StateIn &x ) 235 { 236 return adjust_size_by_resizeability( 237 m_xout, x, typename is_resizeable< state_type >::type() ); 238 } 239 240 template < class StateInOut > extrapolate(size_t k,state_table_type & table,const value_matrix & coeff,StateInOut & xest)241 void extrapolate( size_t k, state_table_type &table, 242 const value_matrix &coeff, StateInOut &xest ) 243 /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf 244 uses the obtained intermediate results to extrapolate to dt->0 245 */ 246 { 247 static const value_type val1 = static_cast< value_type >( 1.0 ); 248 249 for ( int j = k - 1; j > 0; --j ) 250 { 251 stepper_base_type::m_algebra.for_each3( 252 table[j - 1].m_v, table[j].m_v, table[j - 1].m_v, 253 typename operations_type::template scale_sum2< 254 value_type, value_type >( val1 + coeff[k][j], 255 -coeff[k][j] ) ); 256 } 257 stepper_base_type::m_algebra.for_each3( 258 xest, table[0].m_v, xest, 259 typename operations_type::template scale_sum2< 260 value_type, value_type >( val1 + coeff[k][0], -coeff[k][0] ) ); 261 } 262 263 private: 264 midpoint_stepper_type m_midpoint; 265 266 resizer_type m_resizer; 267 resizer_type m_xout_resizer; 268 269 int_vector m_interval_sequence; // stores the successive interval counts 270 value_matrix m_coeff; 271 272 wrapped_state_type m_xout; 273 state_table_type m_table; // sequence of states for extrapolation 274 }; 275 276 /******** DOXYGEN *******/ 277 278 /** 279 * \class extrapolation_stepper 280 * \brief Extrapolation stepper with configurable order, and error estimation. 281 * 282 * The extrapolation stepper is a stepper with error estimation and configurable 283 * order. The order is given as template parameter and needs to be an _odd_ 284 * number. The stepper is based on several executions of the modified midpoint 285 * method and a Richardson extrapolation. This is essentially the same technique 286 * as for bulirsch_stoer, but without the variable order. 287 * 288 * \note The Order parameter has to be an even number greater 2. 289 */ 290 } 291 } 292 } 293 #endif 294