1 /* 2 boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp 3 4 [begin_description] 5 Calculation of the coefficients for the adaptive adams 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_ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED 16 #define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED 17 18 #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp> 19 20 #include <boost/numeric/odeint/util/state_wrapper.hpp> 21 #include <boost/numeric/odeint/util/is_resizeable.hpp> 22 #include <boost/numeric/odeint/util/resizer.hpp> 23 24 #include <boost/numeric/odeint/util/unwrap_reference.hpp> 25 #include <boost/numeric/odeint/util/bind.hpp> 26 27 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> 28 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> 29 30 #include <boost/array.hpp> 31 32 namespace boost { 33 namespace numeric { 34 namespace odeint { 35 namespace detail { 36 37 template< 38 size_t Steps, 39 class Deriv, 40 class Value = double, 41 class Time = double, 42 class Algebra = typename algebra_dispatcher< Deriv >::algebra_type, 43 class Operations = typename operations_dispatcher< Deriv >::operations_type, 44 class Resizer = initially_resizer 45 > 46 class adaptive_adams_coefficients 47 { 48 public: 49 static const size_t steps = Steps; 50 51 typedef unsigned short order_type; 52 static const order_type order_value = steps; 53 54 typedef Value value_type; 55 typedef Deriv deriv_type; 56 typedef Time time_type; 57 58 typedef state_wrapper< deriv_type > wrapped_deriv_type; 59 typedef rotating_buffer< time_type , steps+1 > time_storage_type; 60 61 typedef Algebra algebra_type; 62 typedef Operations operations_type; 63 typedef Resizer resizer_type; 64 65 typedef adaptive_adams_coefficients< Steps , Deriv , Value , Time , Algebra , Operations , Resizer > aac_type; 66 adaptive_adams_coefficients(const algebra_type & algebra=algebra_type ())67 adaptive_adams_coefficients( const algebra_type &algebra = algebra_type()) 68 :m_eo(1), m_steps_init(1), beta(), phi(), m_ns(0), m_time_storage(), 69 m_algebra(algebra), 70 m_phi_resizer() 71 { 72 for (size_t i=0; i<order_value+2; ++i) 73 { 74 c[i] = 1.0/(i+1); 75 c[c_size+i] = 1.0/((i+1)*(i+2)); 76 } 77 78 g[0] = c[0]; 79 g[1] = c[c_size]; 80 81 beta[0][0] = 1; 82 beta[1][0] = 1; 83 84 gs[0] = 1.0; 85 gs[1] = -1.0/2; 86 gs[2] = -1.0/12; 87 gs[3] = -1.0/24; 88 gs[4] = -19.0/720; 89 gs[5] = -3.0/160; 90 gs[6] = -863.0/60480; 91 gs[7] = -275.0/24192; 92 gs[8] = -33953.0/3628800; 93 gs[9] = 35.0/4436; 94 gs[10] = 40.0/5891; 95 gs[11] = 37.0/6250; 96 gs[12] = 25.0/4771; 97 gs[13] = 40.0/8547; 98 }; 99 predict(time_type t,time_type dt)100 void predict(time_type t, time_type dt) 101 { 102 using std::abs; 103 104 m_time_storage[0] = t; 105 106 if (abs(m_time_storage[0] - m_time_storage[1] - dt) > 1e-16 || m_eo >= m_ns) 107 { 108 m_ns = 0; 109 } 110 else if (m_ns < order_value + 2) 111 { 112 m_ns++; 113 } 114 115 for(size_t i=1+m_ns; i<m_eo+1 && i<m_steps_init; ++i) 116 { 117 time_type diff = m_time_storage[0] - m_time_storage[i]; 118 beta[0][i] = beta[0][i-1]*(m_time_storage[0] + dt - m_time_storage[i-1])/diff; 119 } 120 121 for(size_t i=2+m_ns; i<m_eo+2 && i<m_steps_init+1; ++i) 122 { 123 time_type diff = m_time_storage[0] + dt - m_time_storage[i-1]; 124 for(size_t j=0; j<m_eo+1-i+1; ++j) 125 { 126 c[c_size*i+j] = c[c_size*(i-1)+j] - c[c_size*(i-1)+j+1]*dt/diff; 127 } 128 129 g[i] = c[c_size*i]; 130 } 131 }; 132 do_step(const deriv_type & dxdt,const int o=0)133 void do_step(const deriv_type &dxdt, const int o = 0) 134 { 135 m_phi_resizer.adjust_size( dxdt , detail::bind( &aac_type::template resize_phi_impl< deriv_type > , detail::ref( *this ) , detail::_1 ) ); 136 137 phi[o][0].m_v = dxdt; 138 139 for(size_t i=1; i<m_eo+3 && i<m_steps_init+2 && i<order_value+2; ++i) 140 { 141 if (o == 0) 142 { 143 this->m_algebra.for_each3(phi[o][i].m_v, phi[o][i-1].m_v, phi[o+1][i-1].m_v, 144 typename Operations::template scale_sum2<value_type, value_type>(1.0, -beta[o][i-1])); 145 } 146 else 147 { 148 this->m_algebra.for_each2(phi[o][i].m_v, phi[o][i-1].m_v, 149 typename Operations::template scale_sum1<value_type>(1.0)); 150 } 151 } 152 }; 153 confirm()154 void confirm() 155 { 156 beta.rotate(); 157 phi.rotate(); 158 m_time_storage.rotate(); 159 160 if(m_steps_init < order_value+1) 161 { 162 ++m_steps_init; 163 } 164 }; 165 reset()166 void reset() { m_eo = 1; m_steps_init = 1; }; 167 168 size_t m_eo; 169 size_t m_steps_init; 170 171 rotating_buffer<boost::array<value_type, order_value+1>, 2> beta; // beta[0] = beta(n) 172 rotating_buffer<boost::array<wrapped_deriv_type, order_value+2>, 3> phi; // phi[0] = phi(n+1) 173 boost::array<value_type, order_value + 2> g; 174 boost::array<value_type, 14> gs; 175 176 private: 177 template< class StateType > resize_phi_impl(const StateType & x)178 bool resize_phi_impl( const StateType &x ) 179 { 180 bool resized( false ); 181 182 for(size_t i=0; i<(order_value + 2); ++i) 183 { 184 resized |= adjust_size_by_resizeability( phi[0][i], x, typename is_resizeable<deriv_type>::type() ); 185 resized |= adjust_size_by_resizeability( phi[1][i], x, typename is_resizeable<deriv_type>::type() ); 186 resized |= adjust_size_by_resizeability( phi[2][i], x, typename is_resizeable<deriv_type>::type() ); 187 } 188 return resized; 189 }; 190 191 size_t m_ns; 192 193 time_storage_type m_time_storage; 194 static const size_t c_size = order_value + 2; 195 boost::array<value_type, c_size*c_size> c; 196 197 algebra_type m_algebra; 198 199 resizer_type m_phi_resizer; 200 }; 201 202 } // detail 203 } // odeint 204 } // numeric 205 } // boost 206 207 #endif