stepper/detail/adaptive_adams_coefficients.hpp

0001 /*
0002  boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp
0003
0004  [begin_description]
0005  Calculation of the coefficients for the adaptive adams stepper.
0006  [end_description]
0007
0008  Copyright 2017 Valentin Noah Hartmann
0009
0010  Distributed under the Boost Software License, Version 1.0.
0011  (See accompanying file LICENSE_1_0.txt or
0012  copy at http://www.boost.org/LICENSE_1_0.txt)
0013  */
0014
0015 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED
0016 #define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED
0017
0018 #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
0019
0020 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0021 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0022 #include <boost/numeric/odeint/util/resizer.hpp>
0023
0024 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0025 #include <boost/numeric/odeint/util/bind.hpp>
0026
0027 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
0028 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
0029
0030 #include <boost/array.hpp>
0031
0032 namespace boost {
0033 namespace numeric {
0034 namespace odeint {
0035 namespace detail {
0036
0037 template<
0038 size_t Steps,
0039 class Deriv,
0040 class Value = double,
0041 class Time = double,
0042 class Algebra = typename algebra_dispatcher< Deriv >::algebra_type,
0043 class Operations = typename operations_dispatcher< Deriv >::operations_type,
0044 class Resizer = initially_resizer
0045 >
0046 class adaptive_adams_coefficients
0047 {
0048 public:
0049     static const size_t steps = Steps;
0050
0051     typedef unsigned short order_type;
0052     static const order_type order_value = steps;
0053
0054     typedef Value value_type;
0055     typedef Deriv deriv_type;
0056     typedef Time time_type;
0057
0058     typedef state_wrapper< deriv_type > wrapped_deriv_type;
0059     typedef rotating_buffer< time_type , steps+1 > time_storage_type;
0060
0061     typedef Algebra algebra_type;
0062     typedef Operations operations_type;
0063     typedef Resizer resizer_type;
0064
0065     typedef adaptive_adams_coefficients< Steps , Deriv , Value , Time , Algebra , Operations , Resizer > aac_type;
0066
0067     adaptive_adams_coefficients( const algebra_type &algebra = algebra_type())
0068     :m_eo(1), m_steps_init(1), beta(), phi(), m_ns(0), m_time_storage(),
0069     m_algebra(algebra),
0070     m_phi_resizer()
0071     {
0072         for (size_t i=0; i<order_value+2; ++i)
0073         {
0074             c[i] = 1.0/(i+1);
0075             c[c_size+i] = 1.0/((i+1)*(i+2));
0076         }
0077
0078         g[0] = c[0];
0079         g[1] = c[c_size];
0080
0081         beta[0][0] = 1;
0082         beta[1][0] = 1;
0083
0084         gs[0] = 1.0;
0085         gs[1] = -1.0/2;
0086         gs[2] = -1.0/12;
0087         gs[3] = -1.0/24;
0088         gs[4] = -19.0/720;
0089         gs[5] = -3.0/160;
0090         gs[6] = -863.0/60480;
0091         gs[7] = -275.0/24192;
0092         gs[8] = -33953.0/3628800;
0093         gs[9] = 35.0/4436;
0094         gs[10] =  40.0/5891;
0095         gs[11] = 37.0/6250;
0096         gs[12] = 25.0/4771;
0097         gs[13] = 40.0/8547;
0098     };
0099
0100     void predict(time_type t, time_type dt)
0101     {
0102         using std::abs;
0103
0104         m_time_storage[0] = t;
0105
0106         if (abs(m_time_storage[0] - m_time_storage[1] - dt) > 1e-16 || m_eo >= m_ns)
0107         {
0108             m_ns = 0;
0109         }
0110         else if (m_ns < order_value + 2)
0111         {
0112             m_ns++;
0113         }
0114
0115         for(size_t i=1+m_ns; i<m_eo+1 && i<m_steps_init; ++i)
0116         {
0117             time_type diff = m_time_storage[0] - m_time_storage[i];
0118             beta[0][i] = beta[0][i-1]*(m_time_storage[0] + dt - m_time_storage[i-1])/diff;
0119         }
0120
0121         for(size_t i=2+m_ns; i<m_eo+2 && i<m_steps_init+1; ++i)
0122         {
0123             time_type diff = m_time_storage[0] + dt - m_time_storage[i-1];
0124             for(size_t j=0; j<m_eo+1-i+1; ++j)
0125             {
0126                 c[c_size*i+j] = c[c_size*(i-1)+j] - c[c_size*(i-1)+j+1]*dt/diff;
0127             }
0128
0129             g[i] = c[c_size*i];
0130         }
0131     };
0132
0133     void do_step(const deriv_type &dxdt, const int o = 0)
0134     {
0135         m_phi_resizer.adjust_size( dxdt , detail::bind( &aac_type::template resize_phi_impl< deriv_type > , detail::ref( *this ) , detail::_1 ) );
0136
0137         phi[o][0].m_v = dxdt;
0138
0139         for(size_t i=1; i<m_eo+3 && i<m_steps_init+2 && i<order_value+2; ++i)
0140         {
0141             if (o == 0)
0142             {
0143                 this->m_algebra.for_each3(phi[o][i].m_v, phi[o][i-1].m_v, phi[o+1][i-1].m_v,
0144                     typename Operations::template scale_sum2<value_type, value_type>(1.0, -beta[o][i-1]));
0145             }
0146             else
0147             {
0148                 this->m_algebra.for_each2(phi[o][i].m_v, phi[o][i-1].m_v,
0149                     typename Operations::template scale_sum1<value_type>(1.0));
0150             }
0151         }
0152     };
0153
0154     void confirm()
0155     {
0156         beta.rotate();
0157         phi.rotate();
0158         m_time_storage.rotate();
0159
0160         if(m_steps_init < order_value+1)
0161         {
0162             ++m_steps_init;
0163         }
0164     };
0165
0166     void reset() { m_eo = 1; m_steps_init = 1; };
0167
0168     size_t m_eo;
0169     size_t m_steps_init;
0170
0171     rotating_buffer<boost::array<value_type, order_value+1>, 2> beta; // beta[0] = beta(n)
0172     rotating_buffer<boost::array<wrapped_deriv_type, order_value+2>, 3> phi; // phi[0] = phi(n+1)
0173     boost::array<value_type, order_value + 2> g;
0174     boost::array<value_type, 14> gs;
0175
0176 private:
0177     template< class StateType >
0178     bool resize_phi_impl( const StateType &x )
0179     {
0180         bool resized( false );
0181
0182         for(size_t i=0; i<(order_value + 2); ++i)
0183         {
0184             resized |= adjust_size_by_resizeability( phi[0][i], x, typename is_resizeable<deriv_type>::type() );
0185             resized |= adjust_size_by_resizeability( phi[1][i], x, typename is_resizeable<deriv_type>::type() );
0186             resized |= adjust_size_by_resizeability( phi[2][i], x, typename is_resizeable<deriv_type>::type() );
0187         }
0188         return resized;
0189     };
0190
0191     size_t m_ns;
0192
0193     time_storage_type m_time_storage;
0194     static const size_t c_size = order_value + 2;
0195     boost::array<value_type, c_size*c_size> c;
0196
0197     algebra_type m_algebra;
0198
0199     resizer_type m_phi_resizer;
0200 };
0201
0202 } // detail
0203 } // odeint
0204 } // numeric
0205 } // boost
0206
0207 #endif