1 /* 2 boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp 3 4 [begin_description] 5 Implementation of the stepsize controller for the controlled adams bashforth moulton stepper. 6 [end_description] 7 8 Copyright 2017 Valentin Noah Hartmann 9 10 Distributed under the Boost Software License, Version 1.0. 11 (See accompanying file LICENSE_1_0.txt or 12 copy at http://www.boost.org/LICENSE_1_0.txt) 13 */ 14 15 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_HPP_INCLUDED 16 #define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_HPP_INCLUDED 17 18 #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp> 19 #include <boost/numeric/odeint/stepper/detail/pid_step_adjuster_coefficients.hpp> 20 21 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> 22 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> 23 24 #include <math.h> 25 26 namespace boost { 27 namespace numeric { 28 namespace odeint { 29 namespace detail { 30 31 template< 32 class Value = double, 33 class Time = double 34 > 35 struct pid_op 36 { 37 public: 38 typedef Value value_type; 39 typedef Time time_type; 40 41 const double beta1; 42 const double beta2; 43 const double beta3; 44 const double alpha1; 45 const double alpha2; 46 47 const time_type dt1; 48 const time_type dt2; 49 const time_type dt3; 50 51 const size_t m_steps; 52 pid_opboost::numeric::odeint::detail::pid_op53 pid_op(const size_t steps, const double _dt1, const double _dt2, const double _dt3, 54 const double b1 = 1, const double b2 = 0, const double b3 = 0, const double a1 = 0, const double a2 = 0) 55 :beta1(b1), beta2(b2), beta3(b3), alpha1(a1), alpha2(a2), 56 dt1(_dt1), dt2(_dt2), dt3(_dt3), 57 m_steps(steps) 58 {}; 59 60 template<class T1, class T2, class T3, class T4> operator ()boost::numeric::odeint::detail::pid_op61 void operator()(T1 &t1, const T2 &t2, const T3 &t3, const T4 &t4) 62 { 63 using std::abs; 64 65 t1 = adapted_pow(abs(t2), -beta1/(m_steps + 1)) * 66 adapted_pow(abs(t3), -beta2/(m_steps + 1)) * 67 adapted_pow(abs(t4), -beta3/(m_steps + 1)) * 68 adapted_pow(abs(dt1/dt2), -alpha1/(m_steps + 1))* 69 adapted_pow(abs(dt2/dt3), -alpha2/(m_steps + 1)); 70 71 t1 = 1/t1; 72 }; 73 74 template<class T1, class T2> operator ()boost::numeric::odeint::detail::pid_op75 void operator()(T1 &t1, const T2 &t2) 76 { 77 using std::abs; 78 79 t1 = adapted_pow(abs(t2), -beta1/(m_steps + 1)); 80 81 t1 = 1/t1; 82 }; 83 84 private: 85 template<class T> adapted_powboost::numeric::odeint::detail::pid_op86 inline value_type adapted_pow(T base, double exp) 87 { 88 if(exp == 0) 89 { 90 return 1; 91 } 92 else if (exp > 0) 93 { 94 return pow(base, exp); 95 } 96 else 97 { 98 return 1/pow(base, -exp); 99 } 100 }; 101 }; 102 103 template< 104 class State, 105 class Value = double, 106 class Deriv = State, 107 class Time = double, 108 class Algebra = typename algebra_dispatcher< State >::algebra_type, 109 class Operations = typename operations_dispatcher< Deriv >::operations_type, 110 size_t Type = BASIC 111 > 112 struct pid_step_adjuster 113 { 114 public: thresholdboost::numeric::odeint::detail::pid_step_adjuster115 static double threshold() { return 0.9; }; 116 117 typedef State state_type; 118 typedef Value value_type; 119 typedef Deriv deriv_type; 120 typedef Time time_type; 121 122 typedef Algebra algebra_type; 123 typedef Operations operations_type; 124 125 typedef rotating_buffer<state_type, 3> error_storage_type; 126 typedef rotating_buffer<time_type, 3> time_storage_type; 127 typedef pid_step_adjuster_coefficients<Type> coeff_type; 128 pid_step_adjusterboost::numeric::odeint::detail::pid_step_adjuster129 pid_step_adjuster(double abs_tol = 1e-6, double rel_tol = 1e-6, time_type dtmax = 1.0) 130 :m_dtmax(dtmax), m_error_storage(), m_dt_storage(), m_init(0), 131 m_abs_tol(abs_tol), m_rel_tol(rel_tol) 132 {}; 133 adjust_stepsizeboost::numeric::odeint::detail::pid_step_adjuster134 time_type adjust_stepsize(const size_t steps, time_type dt, state_type &err, const state_type &x, const deriv_type &dxdt) 135 { 136 using std::abs; 137 m_algebra.for_each3( err , x , dxdt , 138 typename operations_type::template rel_error< value_type >( m_abs_tol , m_rel_tol , 1.0 , 1.0 * abs(get_unit_value( dt )) ) ); 139 140 m_error_storage[0] = err; 141 m_dt_storage[0] = dt; 142 143 if(m_init >= 2) 144 { 145 m_algebra.for_each4(err, m_error_storage[0], m_error_storage[1], m_error_storage[2], 146 pid_op<>(steps, m_dt_storage[0], m_dt_storage[1], m_dt_storage[2], 147 m_coeff[0], m_coeff[1], m_coeff[2], m_coeff[3], m_coeff[4])); 148 } 149 else 150 { 151 m_algebra.for_each2(err, m_error_storage[0], 152 pid_op<>(steps, m_dt_storage[0], m_dt_storage[1], m_dt_storage[2], 0.7)); 153 } 154 155 value_type ratio = 1 / m_algebra.norm_inf(err); 156 157 value_type kappa = 1.0; 158 ratio = 1.0 + kappa*atan((ratio - 1) / kappa); 159 160 if(ratio*dt >= m_dtmax) 161 { 162 ratio = m_dtmax / dt; 163 } 164 165 if(ratio >= threshold()) 166 { 167 m_error_storage.rotate(); 168 m_dt_storage.rotate(); 169 170 ++m_init; 171 } 172 else 173 { 174 m_init = 0; 175 } 176 177 return dt * static_cast<time_type>(ratio); 178 }; 179 180 private: 181 algebra_type m_algebra; 182 183 time_type m_dtmax; 184 error_storage_type m_error_storage; 185 time_storage_type m_dt_storage; 186 187 size_t m_init; 188 double m_abs_tol; 189 double m_rel_tol; 190 191 coeff_type m_coeff; 192 }; 193 194 } // detail 195 } // odeint 196 } // numeric 197 } // boost 198 199 #endif