diff --git a/.travis.yml b/.travis.yml index 8d97c457..118d3c0f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -19,7 +19,7 @@ addons: env: - CXXSTD="" - - CXXSTD="cxxflags='-std=c++0x'" + - CXXSTD="-std=c++0x" # For now disable gcc on osx as g++4.8 is not yet available matrix: @@ -28,21 +28,19 @@ matrix: compiler: gcc before_install: -# 1.55: http://sourceforge.net/projects/boost/files/boost/1.55.0/boost_1_55_0.tar.bz2/download -# 1.57: http://sourceforge.net/projects/boost/files/boost/1.57.0/boost_1_57_0.tar.bz2/download # 1.58: http://downloads.sourceforge.net/project/boost/boost/1.58.0/boost_1_58_0.tar.bz2\?r\=\&ts\=1435589970\&use_mirror\=garr - - wget http://downloads.sourceforge.net/project/boost/boost/1.58.0/boost_1_58_0.tar.bz2\?r\=\&ts\=1435589970\&use_mirror\=garr -O /tmp/boost.tar.bz2 +# 1.63: https://sourceforge.net/projects/boost/files/boost/1.63.0/boost_1_63_0.tar.bz2/download?use_mirror=superb-dca2 + - wget https://sourceforge.net/projects/boost/files/boost/1.63.0/boost_1_63_0.tar.bz2/download?use_mirror=superb-dca2 -O /tmp/boost.tar.bz2 - tar jxf /tmp/boost.tar.bz2 - - mv boost_1_58_0 $PWD/boost-trunk -# patch the boost build system - not neccessary with 1.58 anymore -# - patch $PWD/boost-trunk/tools/build/v2/build/toolset.jam toolset.jam.patch + - mv boost_1_63_0 $PWD/boost-trunk - export BOOST_ROOT="$PWD/boost-trunk" - cd $BOOST_ROOT - ./bootstrap.sh - cd $TRAVIS_BUILD_DIR - if [ "$TRAVIS_OS_NAME" = "osx" ] && [ "$CC" = "gcc" ]; then export CC=gcc-4.8; fi + - echo $CC - $CC --version script: - - $BOOST_ROOT/b2 toolset=$CC $CXXSTD \ No newline at end of file + - $BOOST_ROOT/b2 toolset=$CC cxxflags='$CXXSTD' -j2 \ No newline at end of file diff --git a/README.md b/README.md index 735afd52..23e4f88b 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,3 @@ -[![Bitdeli Badge](https://d2weczhvl823v0.cloudfront.net/headmyshoulder/odeint-v2/trend.png)](https://bitdeli.com/free "Bitdeli Badge") +[![Build Status](https://travis-ci.org/headmyshoulder/odeint-v2.svg?branch=master)](https://travis-ci.org/headmyshoulder/odeint-v2) odeint is a highly flexible library for solving ordinary differential equations. diff --git a/include/boost/numeric/odeint.hpp b/include/boost/numeric/odeint.hpp index e0100899..642e46fe 100644 --- a/include/boost/numeric/odeint.hpp +++ b/include/boost/numeric/odeint.hpp @@ -51,6 +51,9 @@ #include +#include +#include + #include #include #include diff --git a/include/boost/numeric/odeint/algebra/multi_array_algebra.hpp b/include/boost/numeric/odeint/algebra/multi_array_algebra.hpp index 0bc476e6..f0ae30e4 100644 --- a/include/boost/numeric/odeint/algebra/multi_array_algebra.hpp +++ b/include/boost/numeric/odeint/algebra/multi_array_algebra.hpp @@ -19,6 +19,8 @@ #define BOOST_NUMERIC_ODEINT_ALGEBRA_MULTI_ARRAY_ALGEBRA_HPP_DEFINED +#include + #include #include #include @@ -121,7 +123,7 @@ struct multi_array_algebra { detail::for_each15( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s13.data() , s14.data() , s15.data() , op ); } - + template< typename S > static typename norm_result_type::type norm_inf( const S &s ) { @@ -129,15 +131,11 @@ struct multi_array_algebra } }; -// template< class T , size_t N > -// struct algebra_dispatcher< boost::array< T , N > > -// { -// typedef array_algebra algebra_type; -// }; - - - - +template< class T , size_t N > +struct algebra_dispatcher< boost::multi_array< T , N > > +{ + typedef multi_array_algebra algebra_type; +}; } // namespace odeint diff --git a/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp b/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp index b4ee5c3b..3783cd81 100644 --- a/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp +++ b/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp @@ -25,42 +25,36 @@ // from odeint // (that is, it lets odeint treat the eigen matrices correctly, knowing // how to add, multiply, compute the norm, etc) - namespace Eigen { - template inline const -typename Eigen::CwiseUnaryOp< - typename Eigen::internal::scalar_add_op< - typename Eigen::internal::traits::Scalar>, - const D > +typename Eigen::CwiseBinaryOp< + internal::scalar_sum_op::Scalar>, + typename DenseBase::ConstantReturnType, + const D> operator+(const typename Eigen::MatrixBase &m, const typename Eigen::internal::traits::Scalar &s) { - return Eigen::CwiseUnaryOp< - typename Eigen::internal::scalar_add_op< - typename Eigen::internal::traits::Scalar>, - const D >(m.derived(),Eigen::internal::scalar_add_op< - typename Eigen::internal::traits::Scalar>(s)); + return CwiseBinaryOp< + internal::scalar_sum_op::Scalar>, + typename DenseBase::ConstantReturnType, + const D>(DenseBase::Constant(m.rows(), m.cols(), s), m.derived()); } template inline const -typename Eigen::CwiseUnaryOp< - typename Eigen::internal::scalar_add_op< - typename Eigen::internal::traits::Scalar>, - const D > +typename Eigen::CwiseBinaryOp< + internal::scalar_sum_op::Scalar>, + typename DenseBase::ConstantReturnType, + const D> operator+(const typename Eigen::internal::traits::Scalar &s, -const typename Eigen::MatrixBase &m) { - return Eigen::CwiseUnaryOp< - typename Eigen::internal::scalar_add_op< - typename Eigen::internal::traits::Scalar>, - const D >(m.derived(),Eigen::internal::scalar_add_op< - typename Eigen::internal::traits::Scalar>(s)); + const typename Eigen::MatrixBase &m) { + return CwiseBinaryOp< + internal::scalar_sum_op::Scalar>, + typename DenseBase::ConstantReturnType, + const D>(DenseBase::Constant(m.rows(), m.cols(), s), m.derived()); } - - template inline const typename Eigen::CwiseBinaryOp< @@ -82,14 +76,9 @@ abs( const Eigen::MatrixBase< D > &m ) { return m.cwiseAbs(); } - - } // end Eigen namespace - - - namespace boost { namespace numeric { namespace odeint { diff --git a/include/boost/numeric/odeint/integrate/integrate_n_steps.hpp b/include/boost/numeric/odeint/integrate/integrate_n_steps.hpp index 7f3a49bd..5cc8aa0e 100644 --- a/include/boost/numeric/odeint/integrate/integrate_n_steps.hpp +++ b/include/boost/numeric/odeint/integrate/integrate_n_steps.hpp @@ -23,6 +23,7 @@ #include #include #include +#include namespace boost { namespace numeric { diff --git a/include/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp b/include/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp index 2f7cc4c6..f3edce19 100644 --- a/include/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp +++ b/include/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp @@ -150,6 +150,12 @@ public : } + void reset(void) + { + m_adams_bashforth.reset(); + } + + private: @@ -175,7 +181,7 @@ private: { m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); m_adams_bashforth.do_step( system , in , t , m_x.m_v , dt ); - m_adams_moulton.do_step( system , in , m_x.m_v , t , out , dt , m_adams_bashforth.step_storage() ); + m_adams_moulton.do_step( system , in , m_x.m_v , t+dt , out , dt , m_adams_bashforth.step_storage() ); } else { @@ -293,6 +299,11 @@ private: * \param dt The step size. */ + /** + * \fn adams_bashforth_moulton::reset( void ) + * \brief Resets the internal buffers of the stepper. + */ + } // odeint } // numeric diff --git a/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp b/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp new file mode 100644 index 00000000..0d3d0c12 --- /dev/null +++ b/include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp @@ -0,0 +1,237 @@ +/* + boost/numeric/odeint/stepper/detail/adaptive_adams_bashforth_moulton.hpp + + [begin_description] + Implemetation of an adaptive adams bashforth moulton stepper. + Used as the stepper for the controlled adams bashforth moulton stepper. + [end_description] + + Copyright 2017 Valentin Noah Hartmann + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#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 + +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_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; + + typedef State state_type; + typedef Value value_type; + typedef Deriv deriv_type; + typedef Time time_type; + + typedef state_wrapper< state_type > wrapped_state_type; + typedef state_wrapper< deriv_type > wrapped_deriv_type; + + typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_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; + + 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() + {}; + + 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, m_xerr.m_v); + boost::numeric::odeint::copy( m_xnew.m_v , inOut); + }; + + 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.m_v); + }; + + template< class System > + void do_step(System system, state_type &inOut, time_type t, time_type dt, state_type &xerr) + { + 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); + boost::numeric::odeint::copy( m_xnew.m_v , inOut); + }; + + 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, xerr); + + 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< class ExplicitStepper, class System > + void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt) + { + reset(); + dt = dt/static_cast< time_type >(order_value); + + m_dxdt_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); + + system( inOut , m_dxdt.m_v , t ); + for( size_t i=0 ; i + void initialize(System system, state_type &inOut, time_type &t, time_type dt) + { + reset(); + dt = dt/static_cast< time_type >(order_value); + + for(size_t i=0; ido_step(system, inOut, t, dt); + t += dt; + }; + }; + + template< class System > + void do_step_impl(System system, const state_type & in, time_type t, state_type & out, time_type &dt, state_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 (m_coeff.m_steps_init == 1) + { + system(in, m_dxdt.m_v, t); + m_coeff.do_step(m_dxdt.m_v, 1); + } + + boost::numeric::odeint::copy( in , out ); + 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, m_coeff.phi[0][eO].m_v, + typename Operations::template scale_sum1(dt*(m_coeff.g[eO]))); + }; + + const coeff_type& coeff() const { return m_coeff; }; + coeff_type & coeff() { return m_coeff; }; + + void reset() { m_coeff.reset(); }; + const deriv_type & dxdt() const { return m_dxdt.m_v; }; + +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() ); + }; + template< class StateType > + bool resize_xerr_impl( const StateType &x ) + { + return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable::type() ); + }; + + 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; + wrapped_state_type m_xerr; +}; + +} // odeint +} // numeric +} // boost + +#endif \ No newline at end of file diff --git a/include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp index 4b908333..02c37492 100644 --- a/include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp +++ b/include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp @@ -98,6 +98,7 @@ public: m_interval_sequence( m_k_max+1 ) , m_coeff( m_k_max+1 ) , m_cost( m_k_max+1 ) , + m_facmin_table( m_k_max+1 ) , m_table( m_k_max ) , STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 ) { @@ -112,21 +113,14 @@ public: else m_cost[i] = m_cost[i-1] + m_interval_sequence[i]; m_coeff[i].resize(i); + m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) ); for( size_t k = 0 ; k < i ; ++k ) { const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] ); m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation } - - // crude estimate of optimal order - - m_current_k_opt = 4; - /* no calculation because log10 might not exist for value_type! - const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 ); - m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast( m_k_max-1 ) , logfact )); - */ } - + reset(); } @@ -282,7 +276,7 @@ public: { m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(m_k_max-1) , static_cast(m_current_k_opt)+1 ); new_h = h_opt[k]; - new_h *= m_cost[m_current_k_opt]/m_cost[k]; + new_h *= static_cast(m_cost[m_current_k_opt])/static_cast(m_cost[k]); } else new_h = h_opt[m_current_k_opt]; break; @@ -344,6 +338,12 @@ public: { m_first = true; m_last_step_rejected = false; + // crude estimate of optimal order + m_current_k_opt = 4; + /* no calculation because log10 might not exist for value_type! + const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 ); + m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast( m_k_max-1 ) , logfact )); + */ } @@ -417,7 +417,7 @@ private: BOOST_USING_STD_MAX(); using std::pow; value_type expo( 1.0/(2*k+1) ); - value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo ); + value_type facmin = m_facmin_table[k]; value_type fac; if (error == 0.0) fac=1.0/facmin; @@ -506,6 +506,7 @@ private: int_vector m_interval_sequence; // stores the successive interval counts value_matrix m_coeff; int_vector m_cost; // costs for interval count + value_vector m_facmin_table; // for precomputed facmin to save pow calls state_table_type m_table; // sequence of states for extrapolation diff --git a/include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp b/include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp index d876ca3d..22eb373b 100644 --- a/include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp +++ b/include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp @@ -110,6 +110,7 @@ public: m_interval_sequence( m_k_max+1 ) , m_coeff( m_k_max+1 ) , m_cost( m_k_max+1 ) , + m_facmin_table( m_k_max+1 ) , m_table( m_k_max ) , m_mp_states( m_k_max+1 ) , m_derivs( m_k_max+1 ) , @@ -125,9 +126,13 @@ public: m_interval_sequence[i] = 2 + 4*i; // 2 6 10 14 ... m_derivs[i].resize( m_interval_sequence[i] ); if( i == 0 ) + { m_cost[i] = m_interval_sequence[i]; - else + } else + { m_cost[i] = m_cost[i-1] + m_interval_sequence[i]; + } + m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) ); m_coeff[i].resize(i); for( size_t k = 0 ; k < i ; ++k ) { @@ -386,6 +391,11 @@ public: } +protected: + + time_type m_max_dt; + + private: template< class StateInOut , class StateVector > @@ -429,7 +439,7 @@ private: using std::pow; value_type expo = static_cast(1)/(m_interval_sequence[k-1]); - value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo ); + value_type facmin = m_facmin_table[k]; value_type fac; if (error == 0.0) fac = static_cast(1)/facmin; @@ -662,8 +672,6 @@ private: default_error_checker< value_type, algebra_type , operations_type > m_error_checker; modified_midpoint_dense_out< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint; - time_type m_max_dt; - bool m_control_interpolation; bool m_last_step_rejected; @@ -692,6 +700,7 @@ private: int_vector m_interval_sequence; // stores the successive interval counts value_matrix m_coeff; int_vector m_cost; // costs for interval count + value_vector m_facmin_table; // for precomputed facmin to save pow calls state_vector_type m_table; // sequence of states for extrapolation diff --git a/include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp b/include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp new file mode 100644 index 00000000..85ee4d01 --- /dev/null +++ b/include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp @@ -0,0 +1,322 @@ +/* + boost/numeric/odeint/stepper/detail/controlled_adams_bashforth_moulton.hpp + + [begin_description] + Implemetation of an controlled adams bashforth moulton stepper. + [end_description] + + Copyright 2017 Valentin Noah Hartmann + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) +*/ + +#ifndef BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_STEPPER_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include + +namespace boost { +namespace numeric { +namespace odeint { + +template< +size_t MaxOrder, +class State, +class Value = double, +class Algebra = typename algebra_dispatcher< State >::algebra_type +> +class default_order_adjuster +{ +public: + typedef State state_type; + typedef Value value_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 ) + {}; + + size_t adjust_order(size_t order, size_t init, boost::array &xerr) + { + using std::abs; + + value_type errc = abs(m_algebra.norm_inf(xerr[2].m_v)); + + value_type errm1 = 3*errc; + value_type errm2 = 3*errc; + + if(order > 2) + { + errm2 = abs(m_algebra.norm_inf(xerr[0].m_v)); + } + if(order >= 2) + { + errm1 = abs(m_algebra.norm_inf(xerr[1].m_v)); + } + + size_t o_new = order; + + if(order == 2 && errm1 <= 0.5*errc) + { + o_new = order - 1; + } + else if(order > 2 && errm2 < errc && errm1 < errc) + { + o_new = order - 1; + } + + if(init < order) + { + return order+1; + } + else if(o_new == order - 1) + { + return order-1; + } + else if(order <= MaxOrder) + { + value_type errp = abs(m_algebra.norm_inf(xerr[3].m_v)); + + if(order > 1 && errm1 < errc && errp) + { + return order-1; + } + else if(order < MaxOrder && errp < (0.5-0.25*order/MaxOrder) * errc) + { + return order+1; + } + } + + return order; + }; +private: + algebra_type m_algebra; +}; + +template< +class ErrorStepper, +class StepAdjuster = detail::pid_step_adjuster< typename ErrorStepper::state_type, + typename ErrorStepper::value_type, + typename ErrorStepper::deriv_type, + typename ErrorStepper::time_type, + typename ErrorStepper::algebra_type, + typename ErrorStepper::operations_type, + detail::H211PI + >, +class OrderAdjuster = default_order_adjuster< ErrorStepper::order_value, + typename ErrorStepper::state_type, + typename ErrorStepper::value_type, + typename ErrorStepper::algebra_type + >, +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; + typedef typename stepper_type::time_type time_type; + + typedef typename stepper_type::algebra_type algebra_type; + typedef typename stepper_type::operations_type operations_type; + 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< wrapped_state_type , 4 > error_storage_type; + + typedef typename stepper_type::coeff_type coeff_type; + typedef controlled_adams_bashforth_moulton< ErrorStepper , StepAdjuster , OrderAdjuster , Resizer > controlled_stepper_type; + + controlled_adams_bashforth_moulton(step_adjuster_type step_adjuster = step_adjuster_type()) + :m_stepper(), + m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(), + m_step_adjuster( step_adjuster ), m_order_adjuster() + {}; + + template< class ExplicitStepper, class System > + void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt) + { + m_stepper.initialize(stepper, system, inOut, t, dt); + }; + + template< class System > + void initialize(System system, state_type &inOut, time_type &t, time_type dt) + { + m_stepper.initialize(system, inOut, t, dt); + }; + + template< class ExplicitStepper, class System > + void initialize_controlled(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type &dt) + { + reset(); + coeff_type &coeff = m_stepper.coeff(); + + m_dxdt_resizer.adjust_size( inOut , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); + + controlled_step_result res = fail; + + for( size_t i=0 ; i + controlled_step_result try_step(System system, state_type & inOut, time_type &t, time_type &dt) + { + m_xnew_resizer.adjust_size( inOut , detail::bind( &controlled_stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); + + controlled_step_result res = try_step(system, inOut, t, m_xnew.m_v, dt); + + if(res == success) + { + boost::numeric::odeint::copy( m_xnew.m_v , inOut); + } + + return res; + }; + + template< class System > + 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_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[2].m_v); + + coeff_type &coeff = m_stepper.coeff(); + + time_type dtPrev = dt; + dt = m_step_adjuster.adjust_stepsize(coeff.m_eo, dt, m_xerr[2].m_v, out, m_stepper.dxdt() ); + + if(dt / dtPrev >= step_adjuster_type::threshold()) + { + system(out, m_dxdt.m_v, t+dtPrev); + + coeff.do_step(m_dxdt.m_v); + coeff.confirm(); + + t += dtPrev; + + size_t eo = coeff.m_eo; + + // estimate errors for next step + double factor = 1; + algebra_type m_algebra; + + m_algebra.for_each2(m_xerr[2].m_v, coeff.phi[1][eo].m_v, + typename operations_type::template scale_sum1(factor*dt*(coeff.gs[eo]))); + + if(eo > 1) + { + m_algebra.for_each2(m_xerr[1].m_v, coeff.phi[1][eo-1].m_v, + typename operations_type::template scale_sum1(factor*dt*(coeff.gs[eo-1]))); + } + if(eo > 2) + { + m_algebra.for_each2(m_xerr[0].m_v, coeff.phi[1][eo-2].m_v, + typename operations_type::template scale_sum1(factor*dt*(coeff.gs[eo-2]))); + } + if(eo < order_value && coeff.m_eo < coeff.m_steps_init-1) + { + m_algebra.for_each2(m_xerr[3].m_v, coeff.phi[1][eo+1].m_v, + typename operations_type::template scale_sum1(factor*dt*(coeff.gs[eo+1]))); + } + + // adjust order + coeff.m_eo = m_order_adjuster.adjust_order(coeff.m_eo, coeff.m_steps_init-1, m_xerr); + + return success; + } + else + { + return fail; + } + }; + + void reset() { m_stepper.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_xerr_impl( const StateType &x ) + { + bool resized( false ); + + for(size_t i=0; i::type() ); + } + return resized; + }; + template< class StateType > + bool resize_xnew_impl( const StateType &x ) + { + return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable::type() ); + }; + + stepper_type m_stepper; + + wrapped_deriv_type m_dxdt; + error_storage_type m_xerr; + wrapped_state_type m_xnew; + + resizer_type m_dxdt_resizer; + resizer_type m_xerr_resizer; + resizer_type m_xnew_resizer; + + step_adjuster_type m_step_adjuster; + order_adjuster_type m_order_adjuster; +}; + +} // odeint +} // numeric +} // boost + +#endif \ No newline at end of file diff --git a/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp b/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp index 8ae627fe..be00827f 100644 --- a/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp +++ b/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp @@ -143,7 +143,7 @@ public: error = max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast( pow( static_cast(5.0) , -static_cast(stepper_order) ) ) , error); - time_type dt_old = dt; + // time_type dt_old = dt; unused variable warning //error too small - increase dt and keep the evolution and limit scaling factor to 5.0 dt *= static_cast(9)/static_cast(10) * pow(error, static_cast(-1) / stepper_order); @@ -163,7 +163,7 @@ public: time_type get_max_dt() { return m_max_dt; } -private: +protected: time_type m_max_dt; }; @@ -411,10 +411,11 @@ public: template< class System , class StateIn , class DerivIn , class StateOut > controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt ) { - if( !m_step_adjuster.check_step_size_limit(dt) ) + unwrapped_step_adjuster &step_adjuster = m_step_adjuster; + if( !step_adjuster.check_step_size_limit(dt) ) { // given dt was above step size limit - adjust and return fail; - dt = m_step_adjuster.get_max_dt(); + dt = step_adjuster.get_max_dt(); return fail; } @@ -428,13 +429,13 @@ public: if( max_rel_err > 1.0 ) { // error too big, decrease step size and reject this step - dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order()); + dt = step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order()); return fail; } else { // otherwise, increase step size and accept t += dt; - dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order()); + dt = step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order()); return success; } } @@ -505,6 +506,7 @@ private: stepper_type m_stepper; error_checker_type m_error_checker; step_adjuster_type m_step_adjuster; + typedef typename unwrap_reference< step_adjuster_type >::type unwrapped_step_adjuster; resizer_type m_dxdt_resizer; resizer_type m_xerr_resizer; @@ -584,7 +586,7 @@ public: const step_adjuster_type &step_adjuster = step_adjuster_type() , const stepper_type &stepper = stepper_type() ) - : m_stepper( stepper ) , m_error_checker( error_checker ) , m_step_adjuster(step_adjuster) , + : m_stepper( stepper ) , m_error_checker( error_checker ) , m_step_adjuster(step_adjuster) , m_first_call( true ) { } @@ -751,10 +753,11 @@ public: controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t , StateOut &out , DerivOut &dxdt_out , time_type &dt ) { - if( !m_step_adjuster.check_step_size_limit(dt) ) + unwrapped_step_adjuster &step_adjuster = m_step_adjuster; + if( !step_adjuster.check_step_size_limit(dt) ) { // given dt was above step size limit - adjust and return fail; - dt = m_step_adjuster.get_max_dt(); + dt = step_adjuster.get_max_dt(); return fail; } @@ -770,12 +773,12 @@ public: if( max_rel_err > 1.0 ) { // error too big, decrease step size and reject this step - dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order()); + dt = step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order()); return fail; } // otherwise, increase step size and accept t += dt; - dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order()); + dt = step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order()); return success; } @@ -903,6 +906,7 @@ private: stepper_type m_stepper; error_checker_type m_error_checker; step_adjuster_type m_step_adjuster; + typedef typename unwrap_reference< step_adjuster_type >::type unwrapped_step_adjuster; resizer_type m_dxdt_resizer; resizer_type m_xerr_resizer; diff --git a/include/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp b/include/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp index 94abc5af..83982dc7 100644 --- a/include/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp +++ b/include/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp @@ -126,7 +126,7 @@ public: m_t_old = m_t; m_t += m_dt; toggle_current_state(); - return std::make_pair( m_t_old , m_dt ); + return std::make_pair( m_t_old , m_t ); } /* diff --git a/include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp b/include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp new file mode 100644 index 00000000..c1029a16 --- /dev/null +++ b/include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp @@ -0,0 +1,207 @@ +/* + boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp + + [begin_description] + Calculation of the coefficients for the adaptive adams stepper. + [end_description] + + Copyright 2017 Valentin Noah Hartmann + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#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 + +#include +#include + +#include +#include + +#include + +namespace boost { +namespace numeric { +namespace odeint { +namespace detail { + +template< +size_t Steps, +class Deriv, +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 +> +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< deriv_type > wrapped_deriv_type; + typedef rotating_buffer< time_type , steps+1 > time_storage_type; + + typedef Algebra algebra_type; + typedef Operations operations_type; + typedef Resizer resizer_type; + + typedef adaptive_adams_coefficients< Steps , Deriv , Value , Time , Algebra , Operations , Resizer > aac_type; + + adaptive_adams_coefficients( const algebra_type &algebra = algebra_type()) + :m_eo(1), m_steps_init(1), beta(), phi(), m_ns(0), m_time_storage(), + m_algebra(algebra), + m_phi_resizer() + { + for (size_t i=0; i 1e-16 || m_eo >= m_ns) + { + 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])); + } + else + { + this->m_algebra.for_each2(phi[o][i].m_v, phi[o][i-1].m_v, + typename Operations::template scale_sum1(1.0)); + } + } + }; + + void confirm() + { + beta.rotate(); + phi.rotate(); + m_time_storage.rotate(); + + if(m_steps_init < order_value+1) + { + ++m_steps_init; + } + }; + + void reset() { m_eo = 1; m_steps_init = 1; }; + + size_t m_eo; + size_t m_steps_init; + + rotating_buffer, 2> beta; // beta[0] = beta(n) + rotating_buffer, 3> phi; // phi[0] = phi(n+1) + boost::array g; + boost::array gs; + +private: + template< class StateType > + bool resize_phi_impl( const StateType &x ) + { + bool resized( false ); + + for(size_t i=0; i<(order_value + 2); ++i) + { + 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; + }; + + size_t m_ns; + + time_storage_type m_time_storage; + static const size_t c_size = order_value + 2; + boost::array c; + + algebra_type m_algebra; + + resizer_type m_phi_resizer; +}; + +} // detail +} // odeint +} // numeric +} // boost + +#endif \ No newline at end of file diff --git a/include/boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp b/include/boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp new file mode 100644 index 00000000..7c0f258b --- /dev/null +++ b/include/boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp @@ -0,0 +1,199 @@ +/* + boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp + + [begin_description] + Implementation of the stepsize controller for the controlled adams bashforth moulton stepper. + [end_description] + + Copyright 2017 Valentin Noah Hartmann + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_HPP_INCLUDED + +#include +#include + +#include +#include + +#include + +namespace boost { +namespace numeric { +namespace odeint { +namespace detail { + +template< +class Value = double, +class Time = double +> +struct pid_op +{ +public: + typedef Value value_type; + typedef Time time_type; + + const double beta1; + const double beta2; + const double beta3; + const double alpha1; + const double alpha2; + + const time_type dt1; + const time_type dt2; + const time_type dt3; + + const size_t m_steps; + + pid_op(const size_t steps, const double _dt1, const double _dt2, const double _dt3, + const double b1 = 1, const double b2 = 0, const double b3 = 0, const double a1 = 0, const double a2 = 0) + :beta1(b1), beta2(b2), beta3(b3), alpha1(a1), alpha2(a2), + dt1(_dt1), dt2(_dt2), dt3(_dt3), + m_steps(steps) + {}; + + template + void operator()(T1 &t1, const T2 &t2, const T3 &t3, const T4 &t4) + { + using std::abs; + + t1 = adapted_pow(abs(t2), -beta1/(m_steps + 1)) * + adapted_pow(abs(t3), -beta2/(m_steps + 1)) * + adapted_pow(abs(t4), -beta3/(m_steps + 1)) * + adapted_pow(abs(dt1/dt2), -alpha1/(m_steps + 1))* + adapted_pow(abs(dt2/dt3), -alpha2/(m_steps + 1)); + + t1 = 1/t1; + }; + + template + void operator()(T1 &t1, const T2 &t2) + { + using std::abs; + + t1 = adapted_pow(abs(t2), -beta1/(m_steps + 1)); + + t1 = 1/t1; + }; + +private: + template + inline value_type adapted_pow(T base, double exp) + { + if(exp == 0) + { + return 1; + } + else if (exp > 0) + { + return pow(base, exp); + } + else + { + return 1/pow(base, -exp); + } + }; +}; + +template< +class State, +class Value = double, +class Deriv = State, +class Time = double, +class Algebra = typename algebra_dispatcher< State >::algebra_type, +class Operations = typename operations_dispatcher< Deriv >::operations_type, +size_t Type = BASIC +> +struct pid_step_adjuster +{ +public: + static double threshold() { return 0.9; }; + + 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 rotating_buffer error_storage_type; + typedef rotating_buffer time_storage_type; + typedef pid_step_adjuster_coefficients coeff_type; + + pid_step_adjuster(double abs_tol = 1e-6, double rel_tol = 1e-6, time_type dtmax = 1.0) + :m_dtmax(dtmax), m_error_storage(), m_dt_storage(), m_init(0), + m_abs_tol(abs_tol), m_rel_tol(rel_tol) + {}; + + time_type adjust_stepsize(const size_t steps, time_type dt, state_type &err, const state_type &x, const deriv_type &dxdt) + { + using std::abs; + m_algebra.for_each3( err , x , dxdt , + typename operations_type::template rel_error< value_type >( m_abs_tol , m_rel_tol , 1.0 , 1.0 * abs(get_unit_value( dt )) ) ); + + m_error_storage[0] = err; + m_dt_storage[0] = dt; + + if(m_init >= 2) + { + m_algebra.for_each4(err, m_error_storage[0], m_error_storage[1], m_error_storage[2], + pid_op<>(steps, m_dt_storage[0], m_dt_storage[1], m_dt_storage[2], + m_coeff[0], m_coeff[1], m_coeff[2], m_coeff[3], m_coeff[4])); + } + else + { + m_algebra.for_each2(err, m_error_storage[0], + pid_op<>(steps, m_dt_storage[0], m_dt_storage[1], m_dt_storage[2], 0.7)); + } + + value_type ratio = 1 / m_algebra.norm_inf(err); + + value_type kappa = 1.0; + ratio = 1.0 + kappa*atan((ratio - 1) / kappa); + + if(ratio*dt >= m_dtmax) + { + ratio = m_dtmax / dt; + } + + if(ratio >= threshold()) + { + m_error_storage.rotate(); + m_dt_storage.rotate(); + + ++m_init; + } + else + { + m_init = 0; + } + + return dt * static_cast(ratio); + }; + +private: + algebra_type m_algebra; + + time_type m_dtmax; + error_storage_type m_error_storage; + time_storage_type m_dt_storage; + + size_t m_init; + double m_abs_tol; + double m_rel_tol; + + coeff_type m_coeff; +}; + +} // detail +} // odeint +} // numeric +} // boost + +#endif \ No newline at end of file diff --git a/include/boost/numeric/odeint/stepper/detail/pid_step_adjuster_coefficients.hpp b/include/boost/numeric/odeint/stepper/detail/pid_step_adjuster_coefficients.hpp new file mode 100644 index 00000000..2dee6ed1 --- /dev/null +++ b/include/boost/numeric/odeint/stepper/detail/pid_step_adjuster_coefficients.hpp @@ -0,0 +1,180 @@ +/* + boost/numeric/odeint/stepper/detail/pid_step_adjuster_coefficients.hpp + + [begin_description] + Coefficients for the PID stepsize controller. + [end_description] + + Copyright 2017 Valentin Noah Hartmann + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_COEFFICIENTS_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_COEFFICIENTS_HPP_INCLUDED + +#include + +namespace boost { +namespace numeric { +namespace odeint { +namespace detail { + +enum adjuster_type{ + BASIC, + H0211, + H211b, + H211PI, + H0312, + H312b, + H312PID, + H0321, + H321 +}; + +template +class pid_step_adjuster_coefficients; + +template<> +class pid_step_adjuster_coefficients : public boost::array +{ +public: + pid_step_adjuster_coefficients() + : boost::array() + { + (*this)[0] = 1.0; + (*this)[1] = 0.0; + (*this)[2] = 0.0; + (*this)[3] = 0.0; + (*this)[4] = 0.0; + } +}; + +template<> +class pid_step_adjuster_coefficients : public boost::array +{ +public: + pid_step_adjuster_coefficients() + : boost::array() + { + (*this)[0] = 1.0 / 2.0; + (*this)[1] = 1.0 / 2.0; + (*this)[2] = 0.0; + (*this)[3] = 1.0 / 2.0; + (*this)[4] = 0.0; + } +}; + +template<> +class pid_step_adjuster_coefficients : public boost::array +{ +public: + pid_step_adjuster_coefficients() + : boost::array() + { + (*this)[0] = 1.0 / 5.0; + (*this)[1] = 2.0 / 5.0; + (*this)[2] = 0.0; + (*this)[3] = 1.0 / 5.0; + (*this)[4] = 0.0; + } +}; + +template<> +class pid_step_adjuster_coefficients : public boost::array +{ +public: + pid_step_adjuster_coefficients() + : boost::array() + { + (*this)[0] = 1.0 / 6.0; + (*this)[1] = 2.0 / 6.0; + (*this)[2] = 0.0; + (*this)[3] = 0.0; + (*this)[4] = 0.0; + } +}; + +template<> +class pid_step_adjuster_coefficients : public boost::array +{ +public: + pid_step_adjuster_coefficients() + : boost::array() + { + (*this)[0] = 1.0 / 4.0; + (*this)[1] = 2.0 / 2.0; + (*this)[2] = 1.0 / 4.0; + (*this)[3] = 3.0 / 4.0; + (*this)[4] = 1.0 / 4.0; + } +}; + +template<> +class pid_step_adjuster_coefficients : public boost::array +{ +public: + pid_step_adjuster_coefficients() + : boost::array() + { + (*this)[0] = 1.0 / 6.0; + (*this)[1] = 2.0 / 6.0; + (*this)[2] = 1.0 / 6.0; + (*this)[3] = 3.0 / 6.0; + (*this)[4] = 1.0 / 6.0; + } +}; + +template<> +class pid_step_adjuster_coefficients : public boost::array +{ +public: + pid_step_adjuster_coefficients() + : boost::array() + { + (*this)[0] = 1.0 / 18.0; + (*this)[1] = 2.0 / 9.0; + (*this)[2] = 1.0 / 18.0; + (*this)[3] = 0.0; + (*this)[4] = 0.0; + } +}; + +template<> +class pid_step_adjuster_coefficients : public boost::array +{ +public: + pid_step_adjuster_coefficients() + : boost::array() + { + (*this)[0] = 5.0 / 4.0; + (*this)[1] = 1.0 / 2.0; + (*this)[2] = -3.0 / 4.0; + (*this)[3] = -1.0 / 4.0; + (*this)[4] = -3.0 / 4.0; + } +}; + +template<> +class pid_step_adjuster_coefficients : public boost::array +{ +public: + pid_step_adjuster_coefficients() + : boost::array() + { + (*this)[0] = 1.0 / 3.0; + (*this)[1] = 1.0 / 18.0; + (*this)[2] = -5.0 / 18.0; + (*this)[3] = -5.0 / 16.0; + (*this)[4] = -1.0 / 6.0; + } +}; + +} // detail +} // odeint +} // numeric +} // boost + +#endif \ No newline at end of file diff --git a/include/boost/numeric/odeint/stepper/generation.hpp b/include/boost/numeric/odeint/stepper/generation.hpp index 90ad2b9a..a051d13b 100644 --- a/include/boost/numeric/odeint/stepper/generation.hpp +++ b/include/boost/numeric/odeint/stepper/generation.hpp @@ -29,6 +29,7 @@ #include #include +#include #include diff --git a/include/boost/numeric/odeint/stepper/generation/generation_controlled_adams_bashforth_moulton.hpp b/include/boost/numeric/odeint/stepper/generation/generation_controlled_adams_bashforth_moulton.hpp new file mode 100644 index 00000000..ca95a1f5 --- /dev/null +++ b/include/boost/numeric/odeint/stepper/generation/generation_controlled_adams_bashforth_moulton.hpp @@ -0,0 +1,59 @@ +/* + boost/numeric/odeint/stepper/detail/generation_controlled_adams_bashforth_moulton.hpp + + [begin_description] + Spezialization of the generation functions for creation of the controlled adams bashforth moulton stepper. + [end_description] + + Copyright 2017 Valentin Noah Hartmann + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#ifndef GENERATION_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED +#define GENERATION_CONTROLLED_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED + +#include +#include +#include + +namespace boost { +namespace numeric { +namespace odeint { + +template< size_t Steps, class State , class Value , class Deriv , class Time , class Algebra , class Operations , class Resizer > +struct get_controller< adaptive_adams_bashforth_moulton< Steps, State , Value , Deriv , Time , Algebra , Operations , Resizer > > +{ + typedef adaptive_adams_bashforth_moulton stepper_type; + typedef controlled_adams_bashforth_moulton< stepper_type > type; +}; + +// controller factory for controlled_adams_bashforth_moulton +template< class Stepper > +struct controller_factory< Stepper , controlled_adams_bashforth_moulton< Stepper > > +{ + typedef Stepper stepper_type; + typedef controlled_adams_bashforth_moulton< stepper_type > controller_type; + typedef typename controller_type::step_adjuster_type step_adjuster_type; + typedef typename stepper_type::value_type value_type; + typedef typename stepper_type::value_type time_type; + + controller_type operator()( value_type abs_error , value_type rel_error , const stepper_type &stepper ) + { + return controller_type(step_adjuster_type(abs_error, rel_error)); + } + + controller_type operator()( value_type abs_error , value_type rel_error , + time_type max_dt, const stepper_type &stepper ) + { + return controller_type( step_adjuster_type(abs_error, rel_error, max_dt)); + } +}; + +} +} +} + +#endif \ No newline at end of file diff --git a/include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp b/include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp index 0e5edd32..5ad1600d 100644 --- a/include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp +++ b/include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp @@ -163,6 +163,8 @@ public: if( m_max_dt != static_cast(0) ) { dt = detail::min_abs(m_max_dt, dt_new); + } else { + dt = dt_new; } m_last_rejected = false; return success; @@ -198,7 +200,7 @@ public: -private: +protected: template< class StateIn > bool resize_m_xerr( const StateIn &x ) diff --git a/include/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp b/include/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp index 6695ba6a..d296246a 100644 --- a/include/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp +++ b/include/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp @@ -41,7 +41,8 @@ class rosenbrock4_dense_output public: typedef ControlledStepper controlled_stepper_type; - typedef typename controlled_stepper_type::stepper_type stepper_type; + typedef typename unwrap_reference< controlled_stepper_type >::type unwrapped_controlled_stepper_type; + typedef typename unwrapped_controlled_stepper_type::stepper_type stepper_type; typedef typename stepper_type::value_type value_type; typedef typename stepper_type::state_type state_type; typedef typename stepper_type::wrapped_state_type wrapped_state_type; @@ -54,7 +55,7 @@ public: typedef rosenbrock4_dense_output< ControlledStepper > dense_output_stepper_type; rosenbrock4_dense_output( const controlled_stepper_type &stepper = controlled_stepper_type() ) - : m_stepper( stepper ) , + : m_stepper( stepper ) , m_x1() , m_x2() , m_current_state_x1( true ) , m_t() , m_t_old() , m_dt() @@ -75,16 +76,17 @@ public: template< class System > std::pair< time_type , time_type > do_step( System system ) { + unwrapped_controlled_stepper_type &stepper = m_stepper; failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails controlled_step_result res = fail; m_t_old = m_t; do { - res = m_stepper.try_step( system , get_current_state() , m_t , get_old_state() , m_dt ); + res = stepper.try_step( system , get_current_state() , m_t , get_old_state() , m_dt ); fail_checker(); // check for overflow of failed steps } while( res == fail ); - m_stepper.stepper().prepare_dense_output(); + stepper.stepper().prepare_dense_output(); this->toggle_current_state(); return std::make_pair( m_t_old , m_t ); } @@ -96,20 +98,23 @@ public: template< class StateOut > void calc_state( time_type t , StateOut &x ) { - m_stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t ); + unwrapped_controlled_stepper_type &stepper = m_stepper; + stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t ); } template< class StateOut > void calc_state( time_type t , const StateOut &x ) { - m_stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t ); + unwrapped_controlled_stepper_type &stepper = m_stepper; + stepper.stepper().calc_state( t , x , get_old_state() , m_t_old , get_current_state() , m_t ); } template< class StateType > void adjust_size( const StateType &x ) { - m_stepper.adjust_size( x ); + unwrapped_controlled_stepper_type &stepper = m_stepper; + stepper.adjust_size( x ); resize_impl( x ); } diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 502ee419..68725822 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -48,6 +48,8 @@ test-suite "odeint" [ run adams_bashforth.cpp ] [ run adams_moulton.cpp ] [ run adams_bashforth_moulton.cpp ] + [ run controlled_adams_bashforth_moulton.cpp ] + [ run adaptive_adams_coefficients.cpp ] [ run generic_stepper.cpp ] [ run generic_error_stepper.cpp ] [ run bulirsch_stoer.cpp ] @@ -82,7 +84,6 @@ test-suite "odeint" [ run integrate_overflow.cpp ] [ compile unwrap_boost_reference.cpp ] [ compile unwrap_reference.cpp : -std=c++0x : unwrap_reference_C++11 ] - [ compile-fail unwrap_reference.cpp : -std=c++98 : unwrap_reference_C++98 ] [ compile std_array.cpp : -std=c++0x ] : valgrind diff --git a/test/adams_bashforth.cpp b/test/adams_bashforth.cpp index f7e5e2bf..7574346c 100644 --- a/test/adams_bashforth.cpp +++ b/test/adams_bashforth.cpp @@ -208,8 +208,13 @@ BOOST_AUTO_TEST_CASE( test_auto_initialization ) adams.do_step( lorenz() , x , 0.0 , x , 0.1 ); BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 2 ) ); + adams.reset(); + adams.do_step( lorenz() , x , 0.0 , x , 0.1 ); - BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 2 ) ); + BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 3 ) ); + + adams.do_step( lorenz() , x , 0.0 , x , 0.1 ); + BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 4 ) ); } BOOST_AUTO_TEST_CASE( test_manual_initialization ) diff --git a/test/adams_bashforth_moulton.cpp b/test/adams_bashforth_moulton.cpp index c821ac79..8e7688bb 100644 --- a/test/adams_bashforth_moulton.cpp +++ b/test/adams_bashforth_moulton.cpp @@ -104,8 +104,8 @@ BOOST_AUTO_TEST_CASE( test_instantiation ) s4.do_step( lorenz() , x , t , dt ); s5.do_step( lorenz() , x , t , dt ); s6.do_step( lorenz() , x , t , dt ); -// s7.do_step( lorenz() , x , t , dt ); -// s8.do_step( lorenz() , x , t , dt ); +// s7.do_step( lorenz() , x , t , dt ); +// s8.do_step( lorenz() , x , t , dt ); } diff --git a/test/adaptive_adams_coefficients.cpp b/test/adaptive_adams_coefficients.cpp new file mode 100644 index 00000000..362458f0 --- /dev/null +++ b/test/adaptive_adams_coefficients.cpp @@ -0,0 +1,120 @@ +#include +#ifdef BOOST_MSVC + #pragma warning(disable:4996) +#endif + +#define BOOST_TEST_MODULE odeint_adaptive_adams_coefficients + +#include + +#include + +#include + +#include +#include +#include + +using namespace boost::unit_test; +using namespace boost::numeric::odeint; + +typedef double value_type; + +BOOST_AUTO_TEST_SUITE( adaptive_adams_coefficients_test ) + +typedef boost::mpl::range_c< size_t , 2 , 10 > vector_of_steps; +BOOST_AUTO_TEST_CASE_TEMPLATE( test_step, step_type, vector_of_steps ) +{ + const static size_t steps = step_type::value; + + typedef std::vector deriv_type; + typedef double time_type; + + typedef detail::adaptive_adams_coefficients aac_type; + + std::vector deriv; + deriv.push_back(-1); + + time_type t = 0.0; + time_type dt = 0.1; + + aac_type coeff; + for(size_t i=0; i v(10); + v[0] = 1.0/1.0; + v[1] = 1.0/2.0; + v[2] = 5.0/12.0; + v[3] = 9.0/24.0; + v[4] = 251.0/720.0; + v[5] = 95.0/288.0; + v[6] = 19087.0/60480.0; + v[7] = 5257.0/17280.0; + v[8] = 5311869667636789.0/18014398509481984.0; + + for(size_t i=0; i deriv_type; + typedef double time_type; + + typedef detail::adaptive_adams_coefficients<3, deriv_type, time_type> aac_type; + aac_type c1; + + deriv_type deriv(1); + deriv[0] = 1.0; + + time_type t = 0.0; + time_type dt = 0.01; + + for(size_t i=0; i<3; ++i) + { + c1.predict(t, dt); + c1.do_step(deriv); + c1.confirm(); + + t+= dt; + + if(c1.m_eo < 3) + c1.m_eo ++; + } + + aac_type c2(c1); + BOOST_CHECK_EQUAL(c1.phi[0][0].m_v[0], c2.phi[0][0].m_v[0]); + BOOST_CHECK(&(c1.phi[0][0].m_v) != &(c2.phi[0][0].m_v)); + + aac_type c3; + deriv_type *p1 = &(c3.phi[0][0].m_v); + + c3 = c1; + BOOST_CHECK(p1 == (&(c3.phi[0][0].m_v))); + BOOST_CHECK_EQUAL(c1.phi[0][0].m_v[0], c3.phi[0][0].m_v[0]); +} + +BOOST_AUTO_TEST_SUITE_END() \ No newline at end of file diff --git a/test/controlled_adams_bashforth_moulton.cpp b/test/controlled_adams_bashforth_moulton.cpp new file mode 100644 index 00000000..40d65720 --- /dev/null +++ b/test/controlled_adams_bashforth_moulton.cpp @@ -0,0 +1,56 @@ +#include +#ifdef BOOST_MSVC + #pragma warning(disable:4996) +#endif + +#define BOOST_TEST_MODULE odeint_controlled_adams_bashforth_moulton + +#include + +#include +#include + +using namespace boost::unit_test; +using namespace boost::numeric::odeint; + +struct const_sys +{ + template< class State , class Deriv , class Value > + void operator()( const State &x , Deriv &dxdt , const Value &dt ) const + { + dxdt[0] = 1; + } +}; + +typedef boost::array< double , 1 > state_type; +typedef double value_type; + +BOOST_AUTO_TEST_SUITE( controlled_adams_bashforth_moulton_test ) + +BOOST_AUTO_TEST_CASE( test_instantiation ) +{ + controlled_adams_bashforth_moulton > s1; + controlled_adams_bashforth_moulton > s2; + controlled_adams_bashforth_moulton > s3; + controlled_adams_bashforth_moulton > s4; + controlled_adams_bashforth_moulton > s5; + controlled_adams_bashforth_moulton > s6; + controlled_adams_bashforth_moulton > s7; + controlled_adams_bashforth_moulton > s8; + controlled_adams_bashforth_moulton > s9; + + state_type x = {{ 10.0 }}; + value_type t = 0.0 , dt = 0.01; + + s1.try_step(const_sys(), x, t, dt); + s2.try_step(const_sys(), x, t, dt); + s3.try_step(const_sys(), x, t, dt); + s4.try_step(const_sys(), x, t, dt); + s5.try_step(const_sys(), x, t, dt); + s6.try_step(const_sys(), x, t, dt); + s7.try_step(const_sys(), x, t, dt); + s8.try_step(const_sys(), x, t, dt); + s9.try_step(const_sys(), x, t, dt); +} + +BOOST_AUTO_TEST_SUITE_END() \ No newline at end of file diff --git a/test/integrate.cpp b/test/integrate.cpp index a90e8b37..07d4617f 100644 --- a/test/integrate.cpp +++ b/test/integrate.cpp @@ -17,6 +17,8 @@ #define BOOST_TEST_MODULE odeint_integrate_functions +#include + #include #include #include @@ -54,6 +56,9 @@ #include #include +#include +#include + #include using namespace boost::unit_test; @@ -237,7 +242,12 @@ class stepper_methods : public mpl::vector< controlled_runge_kutta< runge_kutta_dopri5< state_type > > , controlled_runge_kutta< runge_kutta_fehlberg78< state_type > > , bulirsch_stoer< state_type > , - dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > > + dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > >, + adaptive_adams_bashforth_moulton<3, state_type>, + adaptive_adams_bashforth_moulton<5, state_type>, + adaptive_adams_bashforth_moulton<7, state_type>, + controlled_adams_bashforth_moulton >, + controlled_adams_bashforth_moulton > //bulirsch_stoer_dense_out< state_type > > { }; @@ -274,12 +284,16 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_times_test_case , Stepper, stepper_meth BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_n_steps_test_case , Stepper, stepper_methods ) { - perform_integrate_n_steps_test< Stepper > tester; - tester(); - tester( 200 , 0.01 ); - tester( 200 , 0.01 ); - tester( 200 , 0.01 ); - tester( 200 , -0.01 ); + if(!boost::is_same > >::value && + !boost::is_same > >::value) + { + perform_integrate_n_steps_test< Stepper > tester; + tester(); + tester( 200 , 0.01 ); + tester( 200 , 0.01 ); + tester( 200 , 0.01 ); + tester( 200 , -0.01 ); + } } BOOST_AUTO_TEST_SUITE_END() diff --git a/test/numeric/Jamfile.v2 b/test/numeric/Jamfile.v2 index 1915a325..a4643f6a 100644 --- a/test/numeric/Jamfile.v2 +++ b/test/numeric/Jamfile.v2 @@ -28,6 +28,7 @@ test-suite "odeint" [ run rosenbrock.cpp ] [ run adams_bashforth.cpp ] [ run adams_bashforth_moulton.cpp ] + [ run adaptive_adams_bashforth_moulton.cpp ] [ run abm_time_dependent.cpp ] [ run order_quadrature_formula.cpp ] [ run velocity_verlet.cpp ] diff --git a/test/numeric/adaptive_adams_bashforth_moulton.cpp b/test/numeric/adaptive_adams_bashforth_moulton.cpp new file mode 100644 index 00000000..0932dd1d --- /dev/null +++ b/test/numeric/adaptive_adams_bashforth_moulton.cpp @@ -0,0 +1,113 @@ +/* Boost numeric test of the adams-bashforth steppers test file + + Copyright 2013 Karsten Ahnert + Copyright 2013-2015 Mario Mulansky + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) +*/ + +// disable checked iterator warning for msvc +#include +#ifdef BOOST_MSVC + #pragma warning(disable:4996) +#endif + +#define BOOST_TEST_MODULE numeric_adaptive_adams_bashforth_moulton + +#include +#include + +#include + +#include + +#include + +#include + +using namespace boost::unit_test; +using namespace boost::numeric::odeint; +namespace mpl = boost::mpl; + +typedef double value_type; + +typedef boost::array< double , 2 > state_type; +typedef runge_kutta_fehlberg78 initializing_stepper; + +// harmonic oscillator, analytic solution x[0] = sin( t ) +struct osc +{ + void operator()( const state_type &x , state_type &dxdt , const double t ) const + { + dxdt[0] = x[1]; + dxdt[1] = -x[0]; + } +}; + +BOOST_AUTO_TEST_SUITE( numeric_adaptive_adams_bashforth_moulton_test ) + + +/* generic test for all adams bashforth steppers */ +template< class Stepper > +struct perform_adaptive_adams_bashforth_moulton_test +{ + void operator()( void ) + { + Stepper stepper; + initializing_stepper init_stepper; + + const int o = stepper.order()+1; //order of the error is order of approximation + 1 + + const state_type x0 = {{ 0.0 , 1.0 }}; + state_type x1 = x0; + double t = 0.0; + double dt = 0.25; + // initialization, does a number of steps to self-start the stepper with a small stepsize + stepper.initialize( init_stepper, osc() , x1 , t , dt); + double A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] ); + double phi = std::asin(x1[0]/A) - t; + + // now we do the actual step + stepper.do_step( osc() , x1 , t , dt ); + // only examine the error of the adams-bashforth step, not the initialization + const double f = 2.0 * std::abs( A*sin(t+dt+phi) - x1[0] ) / std::pow( dt , o ); // upper bound + + std::cout << o << " , " << f << std::endl; + + /* as long as we have errors above machine precision */ + while( f*std::pow( dt , o ) > 1E-16 ) + { + x1 = x0; + t = 0.0; + stepper.initialize( init_stepper, osc() , x1 , t , dt ); + A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] ); + phi = std::asin(x1[0]/A) - t; + // now we do the actual step + stepper.do_step( osc() , x1 , t , dt ); + stepper.reset(); + // only examine the error of the adams-bashforth step, not the initialization + std::cout << "Testing dt=" << dt << " , " << std::abs( A*sin(t+dt+phi) - x1[0] ) << std::endl; + BOOST_CHECK_LT( std::abs( A*sin(t+dt+phi) - x1[0] ) , f*std::pow( dt , o ) ); + dt *= 0.5; + } + } +}; + +typedef mpl::vector< + adaptive_adams_bashforth_moulton< 2 , state_type > , + adaptive_adams_bashforth_moulton< 3 , state_type > , + adaptive_adams_bashforth_moulton< 4 , state_type > , + adaptive_adams_bashforth_moulton< 5 , state_type > , + adaptive_adams_bashforth_moulton< 6 , state_type > , + adaptive_adams_bashforth_moulton< 7 , state_type > + > adaptive_adams_bashforth_moulton_steppers; + +BOOST_AUTO_TEST_CASE_TEMPLATE( adaptive_adams_bashforth_moulton_test , Stepper, adaptive_adams_bashforth_moulton_steppers ) +{ + perform_adaptive_adams_bashforth_moulton_test< Stepper > tester; + tester(); +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/test/numeric/order_quadrature_formula.cpp b/test/numeric/order_quadrature_formula.cpp index bdf2ae46..923183ff 100644 --- a/test/numeric/order_quadrature_formula.cpp +++ b/test/numeric/order_quadrature_formula.cpp @@ -20,6 +20,8 @@ #include #include +#include "boost/format.hpp" + #include #include diff --git a/test/regression/Jamfile.v2 b/test/regression/Jamfile.v2 index 25bbf0c1..4f64698b 100644 --- a/test/regression/Jamfile.v2 +++ b/test/regression/Jamfile.v2 @@ -26,5 +26,6 @@ test-suite "odeint" [ run regression_147.cpp ] [ compile regression_149.cpp : -std=c++0x ] [ run regression_168.cpp ] + [ run regression_189.cpp ] : valgrind ; diff --git a/test/regression/regression_189.cpp b/test/regression/regression_189.cpp new file mode 100644 index 00000000..54530072 --- /dev/null +++ b/test/regression/regression_189.cpp @@ -0,0 +1,72 @@ +/* + + [begin_description] + Test case for issue 189: + Controlled Rosenbrock stepper fails to increase step size + [end_description] + + Copyright 2016 Karsten Ahnert + Copyright 2016 Mario Mulansky + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#define BOOST_TEST_MODULE odeint_regression_189 + +#include + +#include +#include +#include + +using namespace boost::numeric::odeint; +namespace phoenix = boost::phoenix; + +typedef boost::numeric::ublas::vector< double > vector_type; +typedef boost::numeric::ublas::matrix< double > matrix_type; + +struct stiff_system +{ + void operator()( const vector_type &x , vector_type &dxdt , double /* t */ ) + { + dxdt[ 0 ] = -101.0 * x[ 0 ] - 100.0 * x[ 1 ]; + dxdt[ 1 ] = x[ 0 ]; + } +}; + +struct stiff_system_jacobi +{ + void operator()( const vector_type & /* x */ , matrix_type &J , const double & /* t */ , vector_type &dfdt ) + { + J( 0 , 0 ) = -101.0; + J( 0 , 1 ) = -100.0; + J( 1 , 0 ) = 1.0; + J( 1 , 1 ) = 0.0; + dfdt[0] = 0.0; + dfdt[1] = 0.0; + } +}; + + +BOOST_AUTO_TEST_CASE( regression_189 ) +{ + vector_type x( 2 , 1.0 ); + + size_t num_of_steps = integrate_const( make_dense_output< rosenbrock4< double > >( 1.0e-6 , 1.0e-6 ) , + std::make_pair( stiff_system() , stiff_system_jacobi() ) , + x , 0.0 , 50.0 , 0.01 , + std::cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" ); + // regression: number of steps should be 74 + size_t num_of_steps_expected = 74; + BOOST_CHECK_EQUAL( num_of_steps , num_of_steps_expected ); + + vector_type x2( 2 , 1.0 ); + + size_t num_of_steps2 = integrate_const( make_dense_output< runge_kutta_dopri5< vector_type > >( 1.0e-6 , 1.0e-6 ) , + stiff_system() , x2 , 0.0 , 50.0 , 0.01 , + std::cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" ); + num_of_steps_expected = 1531; + BOOST_CHECK_EQUAL( num_of_steps2 , num_of_steps_expected ); +} diff --git a/test/rosenbrock4.cpp b/test/rosenbrock4.cpp index f0b79d51..0b7bea4b 100644 --- a/test/rosenbrock4.cpp +++ b/test/rosenbrock4.cpp @@ -127,6 +127,39 @@ BOOST_AUTO_TEST_CASE( test_rosenbrock4_dense_output ) stepper.calc_state( 0.5 * ( tr.first + tr.second ) , x ); } +class rosenbrock4_controller_max_dt_adaptable : public rosenbrock4_controller< rosenbrock4< value_type > > +{ + public: + void set_max_dt(value_type max_dt) + { + m_max_dt = max_dt; + } +}; + +BOOST_AUTO_TEST_CASE( test_rosenbrock4_dense_output_ref ) +{ + typedef rosenbrock4_dense_output< boost::reference_wrapper< rosenbrock4_controller_max_dt_adaptable > > stepper_type; + rosenbrock4_controller_max_dt_adaptable c_stepper; + stepper_type stepper( boost::ref( c_stepper ) ); + + typedef stepper_type::state_type state_type; + typedef stepper_type::value_type stepper_value_type; + typedef stepper_type::deriv_type deriv_type; + typedef stepper_type::time_type time_type; + state_type x( 2 ); + x( 0 ) = 0.0 ; x(1) = 1.0; + stepper.initialize( x , 0.0 , 0.1 ); + std::pair< value_type , value_type > tr = stepper.do_step( std::make_pair( sys() , jacobi() ) ); + stepper.calc_state( 0.5 * ( tr.first + tr.second ) , x ); + + // adapt the maximal step size to a small value (smaller than the step size) and make a step + const double max_dt = 1e-8; + c_stepper.set_max_dt(max_dt); + stepper.do_step( std::make_pair( sys() , jacobi() ) ); + // check if the step was done with the requested maximal step size + BOOST_CHECK_CLOSE(max_dt, stepper.current_time_step(), 1e-14); +} + BOOST_AUTO_TEST_CASE( test_rosenbrock4_copy_dense_output ) { typedef rosenbrock4_controller< rosenbrock4< value_type > > controlled_stepper_type; diff --git a/test/runge_kutta_concepts.cpp b/test/runge_kutta_concepts.cpp index a16b04fd..42235dea 100644 --- a/test/runge_kutta_concepts.cpp +++ b/test/runge_kutta_concepts.cpp @@ -88,6 +88,7 @@ struct perform_stepper_test typedef T vector_space_type; void operator()( void ) const { + using std::abs; vector_space_type x; x = 2.0; Stepper stepper; diff --git a/test/runge_kutta_controlled_concepts.cpp b/test/runge_kutta_controlled_concepts.cpp index 2132b092..7f70cc53 100644 --- a/test/runge_kutta_controlled_concepts.cpp +++ b/test/runge_kutta_controlled_concepts.cpp @@ -87,6 +87,7 @@ struct perform_controlled_stepper_test typedef T vector_space_type; void operator()( void ) const { + using std::abs; vector_space_type x; x = 2.0; ControlledStepper controlled_stepper; @@ -134,6 +135,7 @@ struct perform_controlled_stepper_test< ControlledStepper , vector_space_type > { void operator()( void ) const { + using std::abs; vector_space_type x; x = 2.0; ControlledStepper controlled_stepper; diff --git a/test_external/eigen/Jamfile.v2 b/test_external/eigen/Jamfile.v2 index 489f67ce..ad0cc2d2 100644 --- a/test_external/eigen/Jamfile.v2 +++ b/test_external/eigen/Jamfile.v2 @@ -31,6 +31,5 @@ test-suite "odeint" [ run runge_kutta4.cpp ] [ run runge_kutta_dopri5.cpp ] [ run integrate.cpp ] - [ compile-fail fail_integrate.cpp ] : valgrind ; diff --git a/test_external/eigen/fail_integrate.cpp b/test_external/eigen/fail_integrate.cpp deleted file mode 100644 index 4c71b018..00000000 --- a/test_external/eigen/fail_integrate.cpp +++ /dev/null @@ -1,47 +0,0 @@ -/* - [auto_generated] - fail_integrate.cpp - - [begin_description] - tba. - [end_description] - - Copyright 2009-2012 Karsten Ahnert - Copyright 2009-2012 Mario Mulansky - - Distributed under the Boost Software License, Version 1.0. - (See accompanying file LICENSE_1_0.txt or - copy at http://www.boost.org/LICENSE_1_0.txt) - */ - -#include -#ifdef BOOST_MSVC - #pragma warning(disable:4996) -#endif - -#define BOOST_TEST_MODULE odeint_dummy - -#include -#include - -#include - -#include "dummy_odes.hpp" - -using namespace boost::unit_test; -using namespace boost::numeric::odeint; - - -BOOST_AUTO_TEST_SUITE( eigen_fail_integrate ) - -BOOST_AUTO_TEST_CASE( test ) -{ - typedef Eigen::Matrix< double , 1 , 1 > state_type; - state_type x; - x[0] = 10.0; - double t_start = 0.0 , t_end = 1.0 , dt = 0.1; - integrate( constant_system_functor_standard() , x , t_start , t_end , dt ); -} - -BOOST_AUTO_TEST_SUITE_END() -