diff options
Diffstat (limited to 'boost/numeric/odeint/stepper')
16 files changed, 434 insertions, 33 deletions
diff --git a/boost/numeric/odeint/stepper/adams_bashforth.hpp b/boost/numeric/odeint/stepper/adams_bashforth.hpp index 59ed1c9025..5ff1e8358c 100644 --- a/boost/numeric/odeint/stepper/adams_bashforth.hpp +++ b/boost/numeric/odeint/stepper/adams_bashforth.hpp @@ -37,6 +37,7 @@ #include <boost/numeric/odeint/stepper/stepper_categories.hpp> #include <boost/numeric/odeint/stepper/runge_kutta4.hpp> +#include <boost/numeric/odeint/stepper/extrapolation_stepper.hpp> #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp> @@ -44,11 +45,28 @@ #include <boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp> #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp> +#include <boost/mpl/arithmetic.hpp> +#include <boost/mpl/min_max.hpp> +#include <boost/mpl/equal_to.hpp> + +namespace mpl = boost::mpl; + namespace boost { namespace numeric { namespace odeint { + using mpl::int_; + + /* if N >= 4, returns the smallest even number > N, otherwise returns 4 */ + template < int N > + struct order_helper + : mpl::max< typename mpl::eval_if< + mpl::equal_to< mpl::modulus< int_< N >, int_< 2 > >, + int_< 0 > >, + int_< N >, int_< N + 1 > >::type, + int_< 4 > >::type + { }; template< size_t Steps , @@ -59,7 +77,9 @@ class Time = Value , class Algebra = typename algebra_dispatcher< State >::algebra_type , class Operations = typename operations_dispatcher< State >::operations_type , class Resizer = initially_resizer , -class InitializingStepper = runge_kutta4< State , Value , Deriv , Time , Algebra , Operations, Resizer > +class InitializingStepper = extrapolation_stepper< order_helper<Steps>::value, + State, Value, Deriv, Time, + Algebra, Operations, Resizer > > class adams_bashforth : public algebra_stepper_base< Algebra , Operations > { @@ -181,7 +201,8 @@ public : { if( i != 0 ) m_step_storage.rotate(); sys( x , m_step_storage[0].m_v , t ); - stepper.do_step( system , x , m_step_storage[0].m_v , t , dt ); + stepper.do_step_dxdt_impl( system, x, m_step_storage[0].m_v, t, + dt ); t += dt; } m_steps_initialized = steps; @@ -222,7 +243,8 @@ private: { if( m_steps_initialized != 0 ) m_step_storage.rotate(); sys( in , m_step_storage[0].m_v , t ); - m_initializing_stepper.do_step( system , in , m_step_storage[0].m_v , t , out , dt ); + m_initializing_stepper.do_step_dxdt_impl( + system, in, m_step_storage[0].m_v, t, out, dt ); ++m_steps_initialized; } else diff --git a/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp b/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp index 19b0239990..2f7cc4c6fb 100644 --- a/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp +++ b/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp @@ -49,7 +49,8 @@ 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 Resizer = initially_resizer, +class InitializingStepper = runge_kutta4< State , Value , Deriv , Time , Algebra , Operations, Resizer > > class adams_bashforth_moulton { @@ -71,12 +72,13 @@ public : typedef Operations operations_type; typedef Resizer resizer_type; typedef stepper_tag stepper_category; + typedef InitializingStepper initializing_stepper_type; static const size_t steps = Steps; #ifndef DOXYGEN_SKIP - typedef adams_bashforth< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_bashforth_type; + typedef adams_bashforth< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type, initializing_stepper_type > adams_bashforth_type; typedef adams_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_moulton_type; - typedef adams_bashforth_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > stepper_type; + typedef adams_bashforth_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type , initializing_stepper_type> stepper_type; #endif //DOXYGEN_SKIP typedef unsigned short order_type; static const order_type order_value = steps; @@ -158,7 +160,7 @@ private: { m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); m_adams_bashforth.do_step( system , x , t , m_x.m_v , dt ); - m_adams_moulton.do_step( system , x , m_x.m_v , t , x , dt , m_adams_bashforth.step_storage() ); + m_adams_moulton.do_step( system , x , m_x.m_v , t+dt , x , dt , m_adams_bashforth.step_storage() ); } else { diff --git a/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp b/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp index 3ba11a4fe0..08009dc16b 100644 --- a/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp +++ b/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp @@ -152,6 +152,22 @@ public: /* + * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt ) + * + * this version is needed when this stepper is used for initializing + * multistep stepper like adams-bashforth. Hence we provide an explicitely + * named version that is not disabled. Meant for internal use only. + */ + template < class System, class StateInOut, class DerivIn > + void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt, + time_type t, time_type dt ) + { + this->stepper().do_step_impl( system , x , dxdt , t , x , dt ); + } + + + + /* * Version 3 : do_step( sys , in , t , out , dt ) * * this version does not solve the forwarding problem, boost.range can not be used @@ -181,10 +197,21 @@ public: { this->stepper().do_step_impl( system , in , dxdt , t , out , dt ); } - - - - + + /* + * named Version 4: do_step_dxdt_impl( sys , in , dxdt , t , out, dt ) + * + * this version is needed when this stepper is used for initializing + * multistep stepper like adams-bashforth. Hence we provide an explicitely + * named version that is not disabled. Meant for internal use only. + */ + template < class System, class StateIn, class DerivIn, class StateOut > + void do_step_dxdt_impl( System system, const StateIn &in, + const DerivIn &dxdt, time_type t, StateOut &out, + time_type dt ) + { + this->stepper().do_step_impl( system , in , dxdt , t , out , dt ); + } /* * Version 5 :do_step( sys , x , t , dt , xerr ) diff --git a/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp b/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp index a055c7fa5f..b1d751a0ce 100644 --- a/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp +++ b/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp @@ -151,11 +151,27 @@ public: /* + * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt ) + * + * this version is needed when this stepper is used for initializing + * multistep stepper like adams-bashforth. Hence we provide an explicitely + * named version that is not disabled. Meant for internal use only. + */ + template< class System , class StateInOut , class DerivInOut > + void do_step_dxdt_impl( System system , StateInOut &x , DerivInOut &dxdt , time_type t , time_type dt ) + { + m_first_call = true; + this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt ); + } + + /* * version 3 : do_step( sys , in , t , out , dt ) * - * this version does not solve the forwarding problem, boost.range can not be used + * this version does not solve the forwarding problem, boost.range can not + * be used. * - * the disable is needed to avoid ambiguous overloads if state_type = time_type + * the disable is needed to avoid ambiguous overloads if + * state_type = time_type */ template< class System , class StateIn , class StateOut > typename boost::disable_if< boost::is_same< StateIn , time_type > , void >::type @@ -174,12 +190,14 @@ public: * * this version does not solve the forwarding problem, boost.range can not be used */ - template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut > - void do_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type t , - StateOut &out , DerivOut &dxdt_out , time_type dt ) + template< class System, class StateIn, class DerivIn, class StateOut, + class DerivOut > + void do_step( System system, const StateIn &in, const DerivIn &dxdt_in, + time_type t, StateOut &out, DerivOut &dxdt_out, time_type dt ) { m_first_call = true; - this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt ); + this->stepper().do_step_impl( system, in, dxdt_in, t, out, dxdt_out, + dt ); } diff --git a/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp b/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp index 40aab80399..d81c8c7a71 100644 --- a/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp +++ b/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp @@ -139,6 +139,21 @@ public: /* + * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt ) + * + * this version is needed when this stepper is used for initializing + * multistep stepper like adams-bashforth. Hence we provide an explicitely + * named version that is not disabled. Meant for internal use only. + */ + template < class System, class StateInOut, class DerivIn > + void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt, + time_type t, time_type dt ) + { + this->stepper().do_step_impl( system , x , dxdt , t , x , dt ); + } + + + /* * Version 3 : do_step( sys , in , t , out , dt ) * * this version does not solve the forwarding problem, boost.range can not be used @@ -164,6 +179,22 @@ public: this->stepper().do_step_impl( system , in , dxdt , t , out , dt ); } + + /* + * named Version 4: do_step_dxdt_impl( sys , in , dxdt , t , out, dt ) + * + * this version is needed when this stepper is used for initializing + * multistep stepper like adams-bashforth. Hence we provide an explicitely + * named version. Meant for internal use only. + */ + template < class System, class StateIn, class DerivIn, class StateOut > + void do_step_dxdt_impl( System system, const StateIn &in, + const DerivIn &dxdt, time_type t, StateOut &out, + time_type dt ) + { + this->stepper().do_step_impl( system , in , dxdt , t , out , dt ); + } + template< class StateIn > void adjust_size( const StateIn &x ) { diff --git a/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp b/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp index ea3523e5bb..eb09aefc8f 100644 --- a/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp +++ b/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp @@ -182,7 +182,7 @@ private: // stepper for systems with function for dq/dt = f(p) and dp/dt = -f(q) template< class System , class StateIn , class StateOut > - void do_step_impl( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , boost::mpl::true_ ) + void do_step_impl( System system , const StateIn &in , time_type /* t */ , StateOut &out , time_type dt , boost::mpl::true_ ) { typedef typename odeint::unwrap_reference< System >::type system_type; typedef typename odeint::unwrap_reference< typename system_type::first_type >::type coor_deriv_func_type; diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp index b604f86a27..e71008a993 100644 --- a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp +++ b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp @@ -194,7 +194,6 @@ public: static const value_type val1( 1.0 ); - typename odeint::unwrap_reference< System >::type &sys = system; if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) ) { reset(); // system resized -> reset @@ -219,12 +218,12 @@ public: m_midpoint.set_steps( m_interval_sequence[k] ); if( k == 0 ) { - m_midpoint.do_step( sys , in , dxdt , t , out , dt ); + m_midpoint.do_step( system , in , dxdt , t , out , dt ); /* the first step, nothing more to do */ } else { - m_midpoint.do_step( sys , in , dxdt , t , m_table[k-1].m_v , dt ); + m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt ); extrapolate( k , m_table , m_coeff , out ); // get error estimate m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v , @@ -341,7 +340,7 @@ public: resize_m_dxdt( x ); resize_m_xnew( x ); resize_impl( x ); - m_midpoint.adjust_size(); + m_midpoint.adjust_size( x ); } diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp index c352c9605a..0036d1d2a6 100644 --- a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp +++ b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp @@ -155,8 +155,6 @@ public: static const value_type val1( 1.0 ); - typename odeint::unwrap_reference< System >::type &sys = system; - bool reject( true ); time_vector h_opt( m_k_max+1 ); @@ -261,6 +259,7 @@ public: { //calculate dxdt for next step and dense output + typename odeint::unwrap_reference< System >::type &sys = system; sys( out , dxdt_new , t+dt ); //prepare dense output @@ -367,7 +366,7 @@ public: void adjust_size( const StateIn &x ) { resize_impl( x ); - m_midpoint.adjust_size(); + m_midpoint.adjust_size( x ); } diff --git a/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp b/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp index 168bc69202..cc279abbe0 100644 --- a/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp +++ b/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp @@ -117,7 +117,7 @@ struct adams_bashforth_call_algebra< 7 , Algebra , Operations > { //BOOST_ASSERT( false ); // not implemented typedef typename Coefficients::value_type value_type; - Algebra::for_each9( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , + algebra.for_each9( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , typename Operations::template scale_sum8< value_type , Time , Time , Time , Time , Time , Time >( 1.0 , dt * coef[0] , dt * coef[1] , dt * coef[2] , dt * coef[3] , dt * coef[4] , dt * coef[5] , dt * coef[6] ) ); } @@ -132,7 +132,7 @@ struct adams_bashforth_call_algebra< 8 , Algebra , Operations > { //BOOST_ASSERT( false ); // not implemented typedef typename Coefficients::value_type value_type; - Algebra::for_each10( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , steps[7].m_v , + algebra.for_each10( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , steps[7].m_v , typename Operations::template scale_sum9< value_type , Time , Time , Time , Time , Time , Time , Time >( 1.0 , dt * coef[0] , dt * coef[1] , dt * coef[2] , dt * coef[3] , dt * coef[4] , dt * coef[5] , dt * coef[6] , dt * coef[7] ) ); } diff --git a/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp b/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp index f8a1505510..b6f5f2a40b 100644 --- a/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp +++ b/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp @@ -32,7 +32,7 @@ template< class Algebra , class Operations > struct adams_moulton_call_algebra< 1 , Algebra , Operations > { template< class StateIn , class StateOut , class DerivIn , class StepStorage , class Coefficients , class Time > - void operator()( Algebra &algebra , const StateIn &in , StateOut &out , const DerivIn &dxdt , const StepStorage &steps , const Coefficients &coef , Time dt ) const + void operator()( Algebra &algebra , const StateIn &in , StateOut &out , const DerivIn &dxdt , const StepStorage& /* steps */ , const Coefficients &coef , Time dt ) const { typedef typename Coefficients::value_type value_type; algebra.for_each3( out , in , dxdt , typename Operations::template scale_sum2< value_type , Time >( 1.0 , dt * coef[0] ) ); diff --git a/boost/numeric/odeint/stepper/euler.hpp b/boost/numeric/odeint/stepper/euler.hpp index 443f942763..1c7c126b96 100644 --- a/boost/numeric/odeint/stepper/euler.hpp +++ b/boost/numeric/odeint/stepper/euler.hpp @@ -76,7 +76,7 @@ public : { } template< class System , class StateIn , class DerivIn , class StateOut > - void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt ) + void do_step_impl( System /* system */ , const StateIn &in , const DerivIn &dxdt , time_type /* t */ , StateOut &out , time_type dt ) { stepper_base_type::m_algebra.for_each3( out , in , dxdt , typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , dt ) ); @@ -84,7 +84,7 @@ public : } template< class StateOut , class StateIn1 , class StateIn2 > - void calc_state( StateOut &x , time_type t , const StateIn1 &old_state , time_type t_old , const StateIn2 ¤t_state , time_type t_new ) const + void calc_state( StateOut &x , time_type t , const StateIn1 &old_state , time_type t_old , const StateIn2 & /*current_state*/ , time_type /* t_new */ ) const { const time_type delta = t - t_old; stepper_base_type::m_algebra.for_each3( x , old_state , stepper_base_type::m_dxdt.m_v , diff --git a/boost/numeric/odeint/stepper/extrapolation_stepper.hpp b/boost/numeric/odeint/stepper/extrapolation_stepper.hpp new file mode 100644 index 0000000000..3b688e958b --- /dev/null +++ b/boost/numeric/odeint/stepper/extrapolation_stepper.hpp @@ -0,0 +1,293 @@ +/* + [auto_generated] + boost/numeric/odeint/stepper/extrapolation_stepper.hpp + + [begin_description] + extrapolation stepper + [end_description] + + Copyright 2009-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) +*/ + +#ifndef BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED + +#include <iostream> + +#include <algorithm> + +#include <boost/config.hpp> // for min/max guidelines +#include <boost/static_assert.hpp> + +#include <boost/numeric/odeint/util/bind.hpp> +#include <boost/numeric/odeint/util/unwrap_reference.hpp> + +#include <boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp> +#include <boost/numeric/odeint/stepper/modified_midpoint.hpp> +#include <boost/numeric/odeint/stepper/controlled_step_result.hpp> +#include <boost/numeric/odeint/algebra/range_algebra.hpp> +#include <boost/numeric/odeint/algebra/default_operations.hpp> +#include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> +#include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> + +#include <boost/numeric/odeint/util/state_wrapper.hpp> +#include <boost/numeric/odeint/util/is_resizeable.hpp> +#include <boost/numeric/odeint/util/resizer.hpp> +#include <boost/numeric/odeint/util/unit_helper.hpp> +#include <boost/numeric/odeint/util/detail/less_with_sign.hpp> + +namespace boost +{ +namespace numeric +{ +namespace odeint +{ + +template < unsigned short Order, 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 > +#ifndef DOXYGEN_SKIP +class extrapolation_stepper + : public explicit_error_stepper_base< + extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra, + Operations, Resizer >, + Order, Order, Order - 2, State, Value, Deriv, Time, Algebra, + Operations, Resizer > +#else +class extrapolation_stepper : public explicit_error_stepper_base +#endif +{ + + private: + // check for Order being odd + BOOST_STATIC_ASSERT_MSG( + ( ( Order % 2 ) == 0 ) && ( Order > 2 ), + "extrapolation_stepper requires even Order larger than 2" ); + + public: +#ifndef DOXYGEN_SKIP + typedef explicit_error_stepper_base< + extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra, + Operations, Resizer >, + Order, Order, Order - 2, State, Value, Deriv, Time, Algebra, Operations, + Resizer > stepper_base_type; +#else + typedef explicit_error_stepper_base< extrapolation_stepper< ... >, ... > + stepper_base_type; +#endif + + typedef typename stepper_base_type::state_type state_type; + typedef typename stepper_base_type::value_type value_type; + typedef typename stepper_base_type::deriv_type deriv_type; + typedef typename stepper_base_type::time_type time_type; + typedef typename stepper_base_type::algebra_type algebra_type; + typedef typename stepper_base_type::operations_type operations_type; + typedef typename stepper_base_type::resizer_type resizer_type; + +#ifndef DOXYGEN_SKIP + typedef typename stepper_base_type::stepper_type stepper_type; + typedef typename stepper_base_type::wrapped_state_type wrapped_state_type; + typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type; + + typedef std::vector< value_type > value_vector; + typedef std::vector< value_vector > value_matrix; + typedef std::vector< size_t > int_vector; + typedef std::vector< wrapped_state_type > state_table_type; + typedef modified_midpoint< state_type, value_type, deriv_type, time_type, + algebra_type, operations_type, + resizer_type > midpoint_stepper_type; + +#endif // DOXYGEN_SKIP + + typedef unsigned short order_type; + static const order_type order_value = stepper_base_type::order_value; + static const order_type stepper_order_value = + stepper_base_type::stepper_order_value; + static const order_type error_order_value = + stepper_base_type::error_order_value; + + const static size_t m_k_max = ( order_value - 2 ) / 2; + + extrapolation_stepper( const algebra_type &algebra = algebra_type() ) + : stepper_base_type( algebra ), m_interval_sequence( m_k_max + 1 ), + m_coeff( m_k_max + 1 ), m_table( m_k_max ) + { + for ( unsigned short i = 0; i < m_k_max + 1; i++ ) + { + m_interval_sequence[i] = 2 * ( i + 1 ); + m_coeff[i].resize( i ); + 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] = + static_cast< value_type >( 1 ) / + ( r * r - static_cast< value_type >( + 1 ) ); // coefficients for extrapolation + } + } + } + + template < class System, class StateIn, class DerivIn, class StateOut, + class Err > + void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt, + time_type t, StateOut &out, time_type dt, Err &xerr ) + { + // std::cout << "dt: " << dt << std::endl; + // normal step + do_step_impl( system, in, dxdt, t, out, dt ); + + static const value_type val1( 1.0 ); + // additionally, perform the error calculation + stepper_base_type::m_algebra.for_each3( + xerr, out, m_table[0].m_v, + typename operations_type::template scale_sum2< + value_type, value_type >( val1, -val1 ) ); + } + + template < class System, class StateInOut, class DerivIn, class Err > + void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt, + time_type t, time_type dt, Err &xerr ) + { + // normal step + do_step_impl_io( system, inout, dxdt, t, dt ); + + static const value_type val1( 1.0 ); + // additionally, perform the error calculation + stepper_base_type::m_algebra.for_each3( + xerr, inout, m_table[0].m_v, + typename operations_type::template scale_sum2< + value_type, value_type >( val1, -val1 ) ); + } + + template < class System, class StateIn, class DerivIn, class StateOut > + void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt, + time_type t, StateOut &out, time_type dt ) + { + m_resizer.adjust_size( + in, detail::bind( &stepper_type::template resize_impl< StateIn >, + detail::ref( *this ), detail::_1 ) ); + size_t k = 0; + m_midpoint.set_steps( m_interval_sequence[k] ); + m_midpoint.do_step( system, in, dxdt, t, out, dt ); + for ( k = 1; k <= m_k_max; ++k ) + { + m_midpoint.set_steps( m_interval_sequence[k] ); + m_midpoint.do_step( system, in, dxdt, t, m_table[k - 1].m_v, dt ); + extrapolate( k, m_table, m_coeff, out ); + } + } + + template < class System, class StateInOut, class DerivIn > + void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt, + time_type t, time_type dt ) + { + // special care for inout + m_xout_resizer.adjust_size( + inout, + detail::bind( &stepper_type::template resize_m_xout< StateInOut >, + detail::ref( *this ), detail::_1 ) ); + do_step_impl( system, inout, dxdt, t, m_xout.m_v, dt ); + boost::numeric::odeint::copy( m_xout.m_v, inout ); + } + + template < class System, class StateInOut, class DerivIn > + void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt, + time_type t, time_type dt ) + { + do_step_impl_io( system , x , dxdt , t , dt ); + } + + template < class System, class StateIn, class DerivIn, class StateOut > + void do_step_dxdt_impl( System system, const StateIn &in, + const DerivIn &dxdt, time_type t, StateOut &out, + time_type dt ) + { + do_step_impl( system , in , dxdt , t , out , dt ); + } + + + template < class StateIn > void adjust_size( const StateIn &x ) + { + resize_impl( x ); + m_midpoint.adjust_size( x ); + } + + private: + template < class StateIn > bool resize_impl( const StateIn &x ) + { + bool resized( false ); + for ( size_t i = 0; i < m_k_max; ++i ) + resized |= adjust_size_by_resizeability( + m_table[i], x, typename is_resizeable< state_type >::type() ); + return resized; + } + + template < class StateIn > bool resize_m_xout( const StateIn &x ) + { + return adjust_size_by_resizeability( + m_xout, x, typename is_resizeable< state_type >::type() ); + } + + template < class StateInOut > + void extrapolate( size_t k, state_table_type &table, + const value_matrix &coeff, StateInOut &xest ) + /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf + uses the obtained intermediate results to extrapolate to dt->0 + */ + { + static const value_type val1 = static_cast< value_type >( 1.0 ); + + for ( int j = k - 1; j > 0; --j ) + { + stepper_base_type::m_algebra.for_each3( + table[j - 1].m_v, table[j].m_v, table[j - 1].m_v, + typename operations_type::template scale_sum2< + value_type, value_type >( val1 + coeff[k][j], + -coeff[k][j] ) ); + } + stepper_base_type::m_algebra.for_each3( + xest, table[0].m_v, xest, + typename operations_type::template scale_sum2< + value_type, value_type >( val1 + coeff[k][0], -coeff[k][0] ) ); + } + + private: + midpoint_stepper_type m_midpoint; + + resizer_type m_resizer; + resizer_type m_xout_resizer; + + int_vector m_interval_sequence; // stores the successive interval counts + value_matrix m_coeff; + + wrapped_state_type m_xout; + state_table_type m_table; // sequence of states for extrapolation +}; + +/******** DOXYGEN *******/ + +/** + * \class extrapolation_stepper + * \brief Extrapolation stepper with configurable order, and error estimation. + * + * The extrapolation stepper is a stepper with error estimation and configurable + * order. The order is given as template parameter and needs to be an _odd_ + * number. The stepper is based on several executions of the modified midpoint + * method and a Richardson extrapolation. This is essentially the same technique + * as for bulirsch_stoer, but without the variable order. + * + * \note The Order parameter has to be an even number greater 2. + */ +} +} +} +#endif diff --git a/boost/numeric/odeint/stepper/implicit_euler.hpp b/boost/numeric/odeint/stepper/implicit_euler.hpp index f3019889cb..e1c64164bc 100644 --- a/boost/numeric/odeint/stepper/implicit_euler.hpp +++ b/boost/numeric/odeint/stepper/implicit_euler.hpp @@ -144,7 +144,7 @@ private: void solve( state_type &x , matrix_type &m ) { int res = boost::numeric::ublas::lu_factorize( m , m_pm.m_v ); - if( res != 0 ) exit(0); + if( res != 0 ) std::exit(0); boost::numeric::ublas::lu_substitute( m , m_pm.m_v , x ); } diff --git a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp index 4b3824fc08..df4e6c48b6 100644 --- a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp +++ b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp @@ -66,6 +66,7 @@ public: { BOOST_USING_STD_MAX(); using std::abs; + using std::sqrt; const size_t n = x.size(); value_type err = 0.0 , sk = 0.0; diff --git a/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp b/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp index 2b8e70b166..32bda0bd7f 100644 --- a/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp +++ b/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp @@ -121,11 +121,20 @@ public : // dt * m_dxh = k4 sys( m_x_tmp.m_v , m_dxh.m_v , t + dt ); + //x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm ) time_type dt6 = dt / static_cast< value_type >( 6 ); time_type dt3 = dt / static_cast< value_type >( 3 ); stepper_base_type::m_algebra.for_each6( out , in , dxdt , m_dxt.m_v , m_dxm.m_v , m_dxh.m_v , - typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) ); + typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) ); + + // x += dt/6 * m_dxdt + dt/3 * m_dxt ) + // stepper_base_type::m_algebra.for_each4( out , in , dxdt , m_dxt.m_v , + // typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt6 , dt3 ) ); + // // x += dt/3 * m_dxm + dt/6 * m_dxh ) + // stepper_base_type::m_algebra.for_each4( out , out , m_dxm.m_v , m_dxh.m_v , + // typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt3 , dt6 ) ); + } template< class StateType > diff --git a/boost/numeric/odeint/stepper/velocity_verlet.hpp b/boost/numeric/odeint/stepper/velocity_verlet.hpp index 24f3c0e3b6..3a20fc25c9 100644 --- a/boost/numeric/odeint/stepper/velocity_verlet.hpp +++ b/boost/numeric/odeint/stepper/velocity_verlet.hpp @@ -114,7 +114,7 @@ public: algebra_stepper_base_type::m_algebra.for_each4( qout , qin , pin , ain , - typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one_half * dt , one * dt * dt ) ); + typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one * dt , one_half * dt * dt ) ); typename odeint::unwrap_reference< System >::type & sys = system; |