From c5aab453864858338ee780997a49dee7daebb7f7 Mon Sep 17 00:00:00 2001 From: Valentin Hartmann Date: Sun, 2 Jul 2017 15:14:09 +0200 Subject: [PATCH] moving to new formulation --- include/boost/numeric/odeint.hpp | 1 - .../stepper/adaptive_adams_bashforth.hpp | 165 ----------- .../adaptive_adams_bashforth_moulton.hpp | 275 ++++++++++-------- .../odeint/stepper/adaptive_adams_moulton.hpp | 78 ----- .../controlled_adams_bashforth_moulton.hpp | 99 ++++++- .../detail/adaptive_adams_coefficients.hpp | 156 +++++----- .../odeint/stepper/detail/polynomial.hpp | 146 ---------- 7 files changed, 318 insertions(+), 602 deletions(-) delete mode 100644 include/boost/numeric/odeint/stepper/adaptive_adams_bashforth.hpp delete mode 100644 include/boost/numeric/odeint/stepper/adaptive_adams_moulton.hpp delete mode 100644 include/boost/numeric/odeint/stepper/detail/polynomial.hpp diff --git a/include/boost/numeric/odeint.hpp b/include/boost/numeric/odeint.hpp index d577b6a4..642e46fe 100644 --- a/include/boost/numeric/odeint.hpp +++ b/include/boost/numeric/odeint.hpp @@ -51,7 +51,6 @@ #include -#include #include #include diff --git a/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth.hpp b/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth.hpp deleted file mode 100644 index 3b745df1..00000000 --- a/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth.hpp +++ /dev/null @@ -1,165 +0,0 @@ -#ifndef STEPPER_ADAMS_BASHFORTH_HPP_INCLUDED -#define STEPPER_ADAMS_BASHFORTH_HPP_INCLUDED - -#include -#include - -#include -#include -#include - -#include - -#include -#include - -#include -#include -#include - -#include - -namespace boost { -namespace numeric { -namespace odeint { - -template< -size_t Steps, -class State, -class Value = double, -class Deriv = State, -class Time = Value, -class Algebra = typename algebra_dispatcher< State >::algebra_type, -class Operations = typename operations_dispatcher< State >::operations_type , -class Resizer = initially_resizer -> -class adaptive_adams_bashforth: public algebra_stepper_base< Algebra , Operations > -{ -public: - static const size_t steps = Steps; - typedef unsigned short order_type; - static const order_type order_value = steps; - - typedef State state_type; - typedef Value value_type; - typedef Deriv deriv_type; - typedef Time time_type; - typedef Resizer resizer_type; - - typedef Algebra algebra_type; - typedef Operations operations_type; - typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type; - - typedef state_wrapper wrapped_state_type; - typedef state_wrapper wrapped_deriv_type; - typedef stepper_tag stepper_category; - - typedef detail::adaptive_adams_coefficients coeff_type; - - typedef adaptive_adams_bashforth< Steps , State , Value , Deriv , Time , Algebra, Operations, Resizer > stepper_type; - - adaptive_adams_bashforth( const algebra_type &algebra = algebra_type() ) - :algebra_stepper_base_type( algebra ) , - m_dxdt_resizer(), m_xnew_resizer() - {}; - - order_type order() const - { - return order_value; - }; - - template - void do_step(System system, state_type & inOut, time_type t, time_type dt) - { - m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); - - do_step(system, inOut, t, m_xnew.m_v, dt); - boost::numeric::odeint::copy( m_xnew.m_v , inOut); - }; - - template - void do_step(System system, const state_type & in, time_type t, state_type & out, time_type dt) - { - m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); - - system(in, m_dxdt.m_v, t); - m_coeff.step(m_dxdt.m_v, t); - m_coeff.confirm(); - - do_step_impl(m_coeff, in, t, out, dt); - }; - - template - void initialize(System system, state_type &inOut, time_type &t, time_type dt) - { - m_coeff.reset(); - - for(size_t i=0; i0) - { - coeff.poly.add_root(coeff.m_ts[coeff.m_effective_order -1 -i] - coeff.m_ts[0]); - } - - time_type c = coeff.poly.evaluate_integrated(dt); - coeff.m_c[i] = c; - - // predict next state - this->m_algebra.for_each3(out, out, coeff.m_ss[i][coeff.m_effective_order-i-2].m_v, - typename Operations::template scale_sum2(1.0, c)); - } - }; - - const coeff_type& coeff() const - { - return m_coeff; - }; - - coeff_type& coeff() - { - return m_coeff; - }; - - void reset() - { - m_coeff.reset(); - }; - -private: - template< class StateType > - bool resize_dxdt_impl( const StateType &x ) - { - return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable::type() ); - }; - template< class StateType > - bool resize_xnew_impl( const StateType &x ) - { - return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable::type() ); - }; - - resizer_type m_dxdt_resizer; - resizer_type m_xnew_resizer; - - coeff_type m_coeff; - wrapped_deriv_type m_dxdt; - wrapped_state_type m_xnew; -}; - -} -} -} - -#endif diff --git a/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp b/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp index 6d3808c9..afec1b80 100644 --- a/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp +++ b/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp @@ -1,44 +1,46 @@ -#ifndef STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED -#define STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED +#ifndef BOOST_NUMERIC_ODEINT_STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED -#include -#include -#include #include -#include +#include +#include + +#include +#include +#include +#include -#include #include #include #include -#include -#include - #include -#include + +#include +#include + +#include namespace boost { namespace numeric { namespace odeint { - - template< size_t Steps, class State, class Value = double, class Deriv = State, class Time = Value, -class Algebra = typename algebra_dispatcher< State >::algebra_type, +class Algebra = typename algebra_dispatcher< State >::algebra_type , class Operations = typename operations_dispatcher< State >::operations_type , -class Resizer = initially_resizer +class Resizer = initially_resizer > class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , Operations > { public: static const size_t steps = Steps; + typedef unsigned short order_type; static const order_type order_value = steps; @@ -46,132 +48,143 @@ public: typedef Value value_type; typedef Deriv deriv_type; typedef Time time_type; - typedef Resizer resizer_type; - typedef Algebra algebra_type; - typedef Operations operations_type; + typedef state_wrapper< state_type > wrapped_state_type; + typedef state_wrapper< deriv_type > wrapped_deriv_type; + typedef boost::array error_storage_type; + typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type; - - typedef state_wrapper wrapped_state_type; - typedef state_wrapper wrapped_deriv_type; - - typedef detail::adaptive_adams_coefficients coeff_type; - - typedef detail::rotating_buffer error_storage_type; - - typedef adaptive_adams_bashforth aab_type; - typedef adaptive_adams_moulton aam_type; - typedef adaptive_adams_bashforth_moulton< Steps , State , Value , Deriv , Time, Algebra, Operations, Resizer > stepper_type; - + typedef typename algebra_stepper_base_type::algebra_type algebra_type; + typedef typename algebra_stepper_base_type::operations_type operations_type; + typedef Resizer resizer_type; typedef error_stepper_tag stepper_category; - adaptive_adams_bashforth_moulton(const algebra_type &algebra = algebra_type()) - :algebra_stepper_base_type( algebra ), m_aab(algebra), m_aam(algebra), - m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(), - m_coeff() + typedef detail::adaptive_adams_coefficients< Steps, Deriv, Value, Time, Algebra, Operations, Resizer> coeff_type; + typedef adaptive_adams_bashforth_moulton< Steps , State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_type; + + order_type order() const { return order_value; }; + order_type stepper_order() const {return order_value + 1; }; + order_type error_order() const {return order_value; }; + + adaptive_adams_bashforth_moulton( const algebra_type &algebra = algebra_type() ) + :algebra_stepper_base_type( algebra ), m_coeff(), + m_dxdt_resizer(), m_xnew_resizer(), m_xerr_resizer() {}; - order_type order() - { - return order_value; - } - - order_type stepper_order() - { - return order_value + 1; - } - - order_type error_order() - { - return order_value; - } - - template - void do_step(System system, state_type & inOut, time_type t, time_type &dt) - { - m_xerr_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); - - do_step(system, inOut, t, dt, m_xerr.m_v); - }; - - template - void do_step(System system, const state_type & in, time_type t, state_type & out, time_type &dt) - { - m_xerr_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); - - do_step(system, in, t, out, dt, m_xerr.m_v); - }; - - template - void do_step(System system, state_type & inOut, time_type t, time_type &dt, state_type &xerr) + template< class System > + void do_step(System system, state_type &inOut, time_type t, time_type dt ) { m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); - - do_step(system, inOut, t, m_xnew.m_v, dt, xerr); + + do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr[1].m_v); boost::numeric::odeint::copy( m_xnew.m_v , inOut); }; - template - void do_step(System system, const state_type & in, time_type &t, state_type & out, time_type dt, state_type &xerr) - { - do_step_impl(system, m_coeff, in, t, out, dt, xerr); - - system(out, m_dxdt.m_v, t+dt); - m_coeff.step(m_dxdt.m_v, t+dt); - m_coeff.confirm(); + template< class System > + void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt ) + { + do_step(system, in, t, out, dt, m_xerr); }; - template - void do_step_impl(System system, coeff_type coeff, const state_type & in, time_type &t, state_type & out, time_type &dt, state_type &xerr) + template< class System > + void do_step(System system, state_type &inOut, time_type t, time_type dt, state_type &xerr) { - m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); + m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); + + do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr[1].m_v); + boost::numeric::odeint::copy( m_xnew.m_v , inOut); + boost::numeric::odeint::copy( m_xerr[1].m_v , xerr); + }; + + template< class System > + void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt , state_type &xerr) + { + do_step_impl(system, in, t, out, dt, m_xerr); + boost::numeric::odeint::copy( m_xerr[1].m_v , xerr); - if(coeff.m_effective_order == 1) + system(out, m_dxdt.m_v, t+dt); + m_coeff.do_step(m_dxdt.m_v); + m_coeff.confirm(); + + if(m_coeff.m_eo < order_value) + m_coeff.m_eo ++; + }; + + template + void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt) + { + m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); + + for( size_t i=0 ; i(-1.0, 1.0)); }; template void initialize(System system, state_type &inOut, time_type &t, time_type dt) { - m_coeff.reset(); - - for(size_t i=0; ido_step(system, inOut, t, dt/static_cast< Time >(order_value)); + t += dt/static_cast< Time >(order_value); + }; + }; + + template + void do_step_impl(System system, const state_type & in, time_type t, state_type & out, time_type &dt, error_storage_type &xerr) + { + size_t eO = m_coeff.m_eo; + + m_xerr_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); + m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); + + m_coeff.predict(t, dt); + if (eO == 1) + { + system(in, m_dxdt.m_v, t); + m_coeff.do_step(m_dxdt.m_v, 1); } - } - const coeff_type& coeff() const - { - return m_coeff; + out = in; + for(size_t i=0; im_algebra.for_each3(out, out, m_coeff.phi[1][i].m_v, + typename Operations::template scale_sum2(1.0, dt*m_coeff.g[i]*m_coeff.beta[0][i])); + } + + system(out, m_dxdt.m_v, t+dt); + m_coeff.do_step(m_dxdt.m_v); + + this->m_algebra.for_each3(out, out, m_coeff.phi[0][eO].m_v, + typename Operations::template scale_sum2(1.0, dt*m_coeff.g[eO])); + + // error for current order + this->m_algebra.for_each2(xerr[1].m_v, m_coeff.phi[0][eO].m_v, + typename Operations::template scale_sum1(dt*(m_coeff.g[eO]-m_coeff.g[eO-1]))); + + // error for order below/above + if (eO > 1) + { + this->m_algebra.for_each2(xerr[0].m_v, m_coeff.phi[0][eO-1].m_v, + typename Operations::template scale_sum1(dt*(m_coeff.g[eO-1]-m_coeff.g[eO-2]))); + } + + this->m_algebra.for_each2(xerr[2].m_v, m_coeff.phi[0][eO+1].m_v, + typename Operations::template scale_sum1(dt*(m_coeff.g[eO+1]-m_coeff.g[eO]))); }; - coeff_type& coeff() - { - return m_coeff; - }; + const coeff_type& coeff() const {return m_coeff;}; + coeff_type & coeff() {return m_coeff;}; - wrapped_deriv_type m_dxdt; + void reset() { m_coeff.reset(); }; private: template< class StateType > @@ -180,29 +193,35 @@ private: return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable::type() ); }; template< class StateType > - bool resize_xerr_impl( const StateType &x ) - { - return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable::type() ); - }; - template< class StateType > bool resize_xnew_impl( const StateType &x ) { return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable::type() ); }; + template< class StateType > + bool resize_xerr_impl( const StateType &x ) + { + bool resized( false ); - aab_type m_aab; - aam_type m_aam; - - resizer_type m_dxdt_resizer; - resizer_type m_xerr_resizer; - resizer_type m_xnew_resizer; - - wrapped_state_type m_xnew; - wrapped_state_type m_xerr; + for(size_t i=0; i<3; ++i) + { + resized |= adjust_size_by_resizeability( m_xerr[i], x, typename is_resizeable::type() ); + } + return resized; + }; coeff_type m_coeff; + + resizer_type m_dxdt_resizer; + resizer_type m_xnew_resizer; + resizer_type m_xerr_resizer; + + wrapped_deriv_type m_dxdt; + wrapped_state_type m_xnew; + error_storage_type m_xerr; }; -}}} +} // odeint +} // numeric +} // boost -#endif +#endif \ No newline at end of file diff --git a/include/boost/numeric/odeint/stepper/adaptive_adams_moulton.hpp b/include/boost/numeric/odeint/stepper/adaptive_adams_moulton.hpp deleted file mode 100644 index 8623b6a5..00000000 --- a/include/boost/numeric/odeint/stepper/adaptive_adams_moulton.hpp +++ /dev/null @@ -1,78 +0,0 @@ -#ifndef STEPPER_ADAMS_MOULTON_HPP_INCLUDED -#define STEPPER_ADAMS_MOULTON_HPP_INCLUDED - -// helper class for controlled_adams_bashforth_moulton - -#include -#include -#include -#include - -#include -#include -#include - -namespace boost { -namespace numeric { -namespace odeint { - -template< -size_t Steps, -class State, -class Value = double, -class Deriv = State, -class Time = Value, -class Algebra = typename algebra_dispatcher< State >::algebra_type, -class Operations = typename operations_dispatcher< State >::operations_type , -class Resizer = initially_resizer -> -class adaptive_adams_moulton: public algebra_stepper_base< Algebra , Operations > -{ -public: - static const size_t steps = Steps; - typedef unsigned short order_type; - static const order_type order_value = steps + 1; - - typedef State state_type; - typedef Value value_type; - typedef Deriv deriv_type; - typedef Time time_type; - - typedef Algebra algebra_type; - typedef Operations operations_type; - typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type; - - typedef detail::adaptive_adams_coefficients coeff_type; - - typedef adaptive_adams_moulton< Steps , State , Value , Deriv , Time , Resizer > stepper_type; - - adaptive_adams_moulton( const algebra_type &algebra = algebra_type()) - :algebra_stepper_base_type( algebra ) - {}; - - order_type order() const - { - return order_value; - }; - - void do_step(coeff_type & coeff, const state_type & in, - time_type t, state_type & out, time_type &dt) - { - coeff.poly.add_root(0); - out = in; - - // integrating - for(size_t i=0; im_algebra.for_each3(out, out, coeff.m_tss[i][coeff.m_effective_order-i-1].m_v, - typename Operations::template scale_sum2(1.0, c)); - } - }; -}; - -} -} -} - -#endif \ No newline at end of file diff --git a/include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp b/include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp index 6b96150f..fa6c1374 100644 --- a/include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp +++ b/include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp @@ -15,18 +15,59 @@ namespace boost{ namespace numeric{ namespace odeint { +template< +size_t MaxOrder, +class State, +class Algebra = typename algebra_dispatcher< State >::algebra_type +> +class default_order_adjuster +{ +public: + typedef State state_type; + typedef state_wrapper< state_type > wrapped_state_type; + + typedef Algebra algebra_type; + + default_order_adjuster(const algebra_type &algebra = algebra_type()) + : m_algebra(algebra) + {}; + + void adjust_order(size_t &order, const boost::array & xerr) + { + if(order > 1 && m_algebra.norm_inf(xerr[0].m_v) < 0.5 * m_algebra.norm_inf(xerr[1].m_v)) + { + --order; + } + else if (order < MaxOrder && m_algebra.norm_inf(xerr[2].m_v) < 0.5 * m_algebra.norm_inf(xerr[1].m_v)) + { + ++order; + } + }; +private: + algebra_type m_algebra; + +}; + template< class ErrorStepper, class StepAdjuster = detail::pid_step_adjuster, + detail::H312PID + >, +class OrderAdjuster = default_order_adjuster, class Resizer = initially_resizer > class controlled_adams_bashforth_moulton { public: typedef ErrorStepper stepper_type; + + static const typename stepper_type::order_type order_value = stepper_type::order_value; + typedef typename stepper_type::state_type state_type; typedef typename stepper_type::value_type value_type; typedef typename stepper_type::deriv_type deriv_type; @@ -37,19 +78,34 @@ public: typedef Resizer resizer_type; typedef StepAdjuster step_adjuster_type; + typedef OrderAdjuster order_adjuster_type; typedef controlled_stepper_tag stepper_category; typedef typename stepper_type::wrapped_state_type wrapped_state_type; typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type; + typedef boost::array error_storage_type; - typedef controlled_adams_bashforth_moulton controlled_stepper_type; + typedef typename stepper_type::coeff_type coeff_type; + typedef controlled_adams_bashforth_moulton controlled_stepper_type; controlled_adams_bashforth_moulton(step_adjuster_type step_adjuster = step_adjuster_type()) - :m_stepper(), m_coeff(m_stepper.coeff()), + :m_stepper(), m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(), - m_step_adjuster(step_adjuster) + m_step_adjuster(step_adjuster), m_order_adjuster() {}; + template + void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt) + { + m_stepper.initialize(stepper, system, inOut, t, dt); + }; + + template + void initialize(System system, state_type &inOut, time_type &t, time_type dt) + { + m_stepper.initialize(system, inOut, t, dt); + }; + template controlled_step_result try_step(System system, state_type & inOut, time_type &t, time_type &dt) { @@ -67,18 +123,23 @@ public: controlled_step_result try_step(System system, const state_type & in, time_type &t, state_type & out, time_type &dt) { m_xerr_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); - m_stepper.do_step_impl(system, m_coeff, in, t, out, dt, m_xerr.m_v); + m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); + + m_stepper.do_step_impl(system, in, t, out, dt, m_xerr); + + coeff_type &coeff = m_stepper.coeff(); + + size_t prevOrder = coeff.m_eo; + m_order_adjuster.adjust_order(coeff.m_eo, m_xerr); time_type dtPrev = dt; - dt = m_step_adjuster.adjust_stepsize(dt, m_xerr.m_v); + dt = m_step_adjuster.adjust_stepsize(dt, m_xerr[1 + coeff.m_eo - prevOrder].m_v); - if(dt / dtPrev >= 0.9) + if(dt / dtPrev >= step_adjuster_type::threshold()) { - m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); - - system(out, m_dxdt.m_v, t+dtPrev); - m_coeff.step(m_dxdt.m_v, t+dtPrev); - m_coeff.confirm(); + system(out, m_dxdt.m_v, t+dt); + coeff.do_step(m_dxdt.m_v); + coeff.confirm(); t += dtPrev; return success; @@ -89,6 +150,8 @@ public: } }; + void reset() { m_stepper.reset(); }; + private: template< class StateType > bool resize_dxdt_impl( const StateType &x ) @@ -98,7 +161,13 @@ private: template< class StateType > bool resize_xerr_impl( const StateType &x ) { - return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable::type() ); + bool resized( false ); + + for(size_t i=0; i<3; ++i) + { + resized |= adjust_size_by_resizeability( m_xerr[i], x, typename is_resizeable::type() ); + } + return resized; }; template< class StateType > bool resize_xnew_impl( const StateType &x ) @@ -107,10 +176,9 @@ private: }; stepper_type m_stepper; - typename stepper_type::coeff_type &m_coeff; wrapped_deriv_type m_dxdt; - wrapped_state_type m_xerr; + error_storage_type m_xerr; wrapped_state_type m_xnew; resizer_type m_dxdt_resizer; @@ -118,6 +186,7 @@ private: resizer_type m_xnew_resizer; step_adjuster_type m_step_adjuster; + order_adjuster_type m_order_adjuster; }; } diff --git a/include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp b/include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp index aad33cfa..f7a32787 100644 --- a/include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp +++ b/include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp @@ -1,8 +1,7 @@ -#ifndef ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED -#define ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED +#ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED #include -#include #include #include @@ -12,8 +11,6 @@ #include #include -#include - #include #include @@ -25,112 +22,133 @@ namespace detail { template< size_t Steps, class Deriv, -class Time, +class Value = double, +class Time = double, class Algebra = typename algebra_dispatcher< Deriv >::algebra_type, class Operations = typename operations_dispatcher< Deriv >::operations_type , class Resizer = initially_resizer > -struct adaptive_adams_coefficients +class adaptive_adams_coefficients { public: static const size_t steps = Steps; + typedef unsigned short order_type; + static const order_type order_value = steps; + + typedef Value value_type; typedef Deriv deriv_type; typedef Time time_type; + typedef state_wrapper wrapped_deriv_type; - typedef detail::rotating_buffer step_storage_type; // +1 for moulton - typedef detail::rotating_buffer time_storage_type; + typedef rotating_buffer step_storage_type; // +1 for moulton + typedef rotating_buffer time_storage_type; typedef Algebra algebra_type; typedef Operations operations_type; - typedef Resizer resizer_type; - typedef adaptive_adams_coefficients aac_type; - typedef detail::Polynomial poly_type; + typedef adaptive_adams_coefficients aac_type; - adaptive_adams_coefficients(const algebra_type &algebra = algebra_type() ) - :poly(), m_effective_order(1), m_resizer(), m_algebra(algebra) - {}; - - void step(const deriv_type &deriv, const time_type &t) + adaptive_adams_coefficients( const algebra_type &algebra = algebra_type()) + :m_eo(1), beta(), phi(), m_ns(0), m_time_storage(), + m_algebra(algebra), + m_phi_resizer() { - m_resizer.adjust_size( deriv , detail::bind( &aac_type::template resize_tss_impl< deriv_type > , detail::ref( *this ) , detail::_1 ) ); + for (size_t i=0; i 1e-16 || m_eo >= m_ns) { - time_type dt = t - m_ts[i-1]; - m_algebra.for_each3(m_tss[i][0].m_v, m_tss[i-1][0].m_v, m_ss[i-1][0].m_v, typename Operations::template scale_sum2(1/dt, -1/dt)); + m_ns = 0; + } + else if (m_ns < order_value + 2) + { + m_ns++; + } + + for(size_t i=1+m_ns; i , detail::ref( *this ) , detail::_1 ) ); + + phi[o][0].m_v = dxdt; + + for(size_t i=1; im_algebra.for_each3(phi[o][i].m_v, phi[o][i-1].m_v, phi[o+1][i-1].m_v, + typename Operations::template scale_sum2(1.0, -beta[o][i-1])); + } + }; + void confirm() { - for(size_t i=0; i m_ss; - time_storage_type m_ts; - - boost::array m_tss; - time_storage_type m_tts; - - size_t m_effective_order; + rotating_buffer, 2> beta; // beta[0] = beta(n) + rotating_buffer, 3> phi; // phi[0] = phi(n+1) + boost::array g; private: template< class StateType > - bool resize_tss_impl( const StateType &x ) + bool resize_phi_impl( const StateType &x ) { + bool resized( false ); - for( size_t i=0 ; i::type() ); - - m_tts[i] = 0; + resized |= adjust_size_by_resizeability( phi[0][i], x, typename is_resizeable::type() ); + resized |= adjust_size_by_resizeability( phi[1][i], x, typename is_resizeable::type() ); + resized |= adjust_size_by_resizeability( phi[2][i], x, typename is_resizeable::type() ); } return resized; }; - resizer_type m_resizer; + size_t m_ns; + + time_storage_type m_time_storage; + boost::array, order_value + 2> c; + algebra_type m_algebra; + + resizer_type m_phi_resizer; }; } @@ -138,4 +156,4 @@ private: } } -#endif +#endif \ No newline at end of file diff --git a/include/boost/numeric/odeint/stepper/detail/polynomial.hpp b/include/boost/numeric/odeint/stepper/detail/polynomial.hpp deleted file mode 100644 index 6ab80e32..00000000 --- a/include/boost/numeric/odeint/stepper/detail/polynomial.hpp +++ /dev/null @@ -1,146 +0,0 @@ -#ifndef POLYNOMIAL_HPP_INCLUDED -#define POLYNOMIAL_HPP_INCLUDED - -#include -#include - -namespace boost { -namespace numeric { -namespace odeint { -namespace detail{ - -template< -size_t order, -class Time -> -class Polynomial -{ -public: - typedef Time time_type; - typedef Polynomial poly_type; - - Polynomial() - { - for(size_t i = 0; i coeff) - { - m_coeff = coeff; - }; - - Polynomial integrate() - { - boost::array coeff_int; - coeff_int[order] = 0; - - for(size_t i=0; i polyInt(coeff_int); - return polyInt; - }; - - time_type evaluate(const time_type &t) - { - // fma: x*y+z - m_res = m_coeff[0]; - for(size_t i=1; i= order)?0:m_coeff[i]/(order-i))); - m_res = ((i >= order)?0:m_coeff[i]/(order-i)) + m_res * t; - } - - return m_res; - } - - void add_root(const time_type &root) - { - for(size_t j=0; j from_roots(time_type * roots, unsigned short num_roots) - { - time_type coeff[order] = {0}; - coeff[order-1] = 1; - - for(size_t i=0; i poly(coeff); - return poly; - }; - - // first element is highest order - boost::array m_coeff; -private: - time_type m_res; -}; - -} -} -} -} - -#endif \ No newline at end of file