boost_mirror/odeint
1
0
Fork 0
mirror of https://github.com/boostorg/odeint.git synced 2025-05-11 21:44:03 +00:00
Code Issues Packages Projects Releases Wiki Activity

changed tabs to spaces in all files

Browse Source
This commit is contained in:
Valentin Hartmann 2017-06-22 17:44:49 +02:00
parent 079f3ffa4d
commit d0dad9a53a
12 changed files with 755 additions and 780 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

172
include/boost/numeric/odeint/stepper/adaptive_adams_bashforth.hpp
View File

@ -35,114 +35,114 @@ class Resizer = initially_resizer
> >
class adaptive_adams_bashforth: public algebra_stepper_base< Algebra , Operations > class adaptive_adams_bashforth: public algebra_stepper_base< Algebra , Operations >
{ {
public: public:
static const size_t steps = Steps; static const size_t steps = Steps;
typedef unsigned short order_type; typedef unsigned short order_type;
static const order_type order_value = steps; static const order_type order_value = steps;
typedef State state_type; typedef State state_type;
typedef Value value_type; typedef Value value_type;
typedef Deriv deriv_type; typedef Deriv deriv_type;
typedef Time time_type; typedef Time time_type;
typedef Resizer resizer_type; typedef Resizer resizer_type;
typedef Algebra algebra_type; typedef Algebra algebra_type;
typedef Operations operations_type; typedef Operations operations_type;
typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type; typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
typedef state_wrapper<state_type> wrapped_state_type; typedef state_wrapper<state_type> wrapped_state_type;
typedef state_wrapper<deriv_type> wrapped_deriv_type; typedef state_wrapper<deriv_type> wrapped_deriv_type;
typedef stepper_tag stepper_category; typedef stepper_tag stepper_category;
typedef detail::adaptive_adams_coefficients<Steps, deriv_type, time_type, algebra_type, operations_type> coeff_type; typedef detail::adaptive_adams_coefficients<Steps, deriv_type, time_type, algebra_type, operations_type> coeff_type;
typedef adaptive_adams_bashforth< Steps , State , Value , Deriv , Time , Algebra, Operations, Resizer > stepper_type; typedef adaptive_adams_bashforth< Steps , State , Value , Deriv , Time , Algebra, Operations, Resizer > stepper_type;
adaptive_adams_bashforth( const algebra_type &algebra = algebra_type() ) adaptive_adams_bashforth( const algebra_type &algebra = algebra_type() )
:algebra_stepper_base_type( algebra ) , :algebra_stepper_base_type( algebra ) ,
m_dxdt_resizer(), m_xnew_resizer() m_dxdt_resizer(), m_xnew_resizer()
{}; {};
order_type order() const order_type order() const
{ {
return order_value; return order_value;
}; };
template<class System> template<class System>
void do_step(System system, state_type & inOut, time_type t, time_type dt) 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 ) ); 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); do_step(system, inOut, t, m_xnew.m_v, dt);
boost::numeric::odeint::copy( m_xnew.m_v , inOut); boost::numeric::odeint::copy( m_xnew.m_v , inOut);
}; };
template<class System> template<class System>
void do_step(System system, const state_type & in, time_type t, state_type & out, time_type dt) void do_step(System system, const state_type & in, time_type t, state_type & out, time_type dt)
{ {
m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
system(in, m_dxdt.m_v, t); system(in, m_dxdt.m_v, t);
m_coeff.step(m_dxdt.m_v, t); m_coeff.step(m_dxdt.m_v, t);
m_coeff.confirm(); m_coeff.confirm();
do_step_impl(m_coeff, in, t, out, dt); do_step_impl(m_coeff, in, t, out, dt);
}; };
void do_step_impl(coeff_type & coeff, const state_type & in, time_type t, state_type & out, time_type dt) void do_step_impl(coeff_type & coeff, const state_type & in, time_type t, state_type & out, time_type dt)
{ {
coeff.poly.reset(); coeff.poly.reset();
out = in; out = in;
// integrating // integrating
for(size_t i=0; i<coeff.m_effective_order-1; ++i) for(size_t i=0; i<coeff.m_effective_order-1; ++i)
{ {
if(i>0) if(i>0)
{ {
coeff.poly.add_root(coeff.m_ts[coeff.m_effective_order -1 -i] - coeff.m_ts[0]); coeff.poly.add_root(coeff.m_ts[coeff.m_effective_order -1 -i] - coeff.m_ts[0]);
} }
time_type c = coeff.poly.evaluate_integrated(dt); time_type c = coeff.poly.evaluate_integrated(dt);
coeff.m_c[i] = c; coeff.m_c[i] = c;
// predict next state // predict next state
this->m_algebra.for_each3(out, out, coeff.m_ss[i][coeff.m_effective_order-i-2].m_v, typename Operations::template scale_sum2<double, double>(1.0, c)); this->m_algebra.for_each3(out, out, coeff.m_ss[i][coeff.m_effective_order-i-2].m_v, typename Operations::template scale_sum2<double, double>(1.0, c));
} }
}; };
const coeff_type& coeff() const const coeff_type& coeff() const
{ {
return m_coeff; return m_coeff;
}; };
coeff_type& coeff() coeff_type& coeff()
{ {
return m_coeff; return m_coeff;
}; };
void reset() void reset()
{ {
m_coeff.reset(); m_coeff.reset();
}; };
private: private:
template< class StateType > template< class StateType >
bool resize_dxdt_impl( const StateType &x ) bool resize_dxdt_impl( const StateType &x )
{ {
return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() ); return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() );
}; };
template< class StateType > template< class StateType >
bool resize_xnew_impl( const StateType &x ) bool resize_xnew_impl( const StateType &x )
{ {
return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() ); return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() );
}; };
resizer_type m_dxdt_resizer; resizer_type m_dxdt_resizer;
resizer_type m_xnew_resizer; resizer_type m_xnew_resizer;
coeff_type m_coeff; coeff_type m_coeff;
wrapped_deriv_type m_dxdt; wrapped_deriv_type m_dxdt;
wrapped_state_type m_xnew; wrapped_state_type m_xnew;
}; };
} }

238
include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp
View File

@ -37,158 +37,158 @@ class Resizer = initially_resizer
> >
class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , Operations > class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , Operations >
{ {
public: public:
static const size_t steps = Steps; static const size_t steps = Steps;
typedef unsigned short order_type; typedef unsigned short order_type;
static const order_type order_value = steps; static const order_type order_value = steps;
typedef State state_type; typedef State state_type;
typedef Value value_type; typedef Value value_type;
typedef Deriv deriv_type; typedef Deriv deriv_type;
typedef Time time_type; typedef Time time_type;
typedef Resizer resizer_type; typedef Resizer resizer_type;
typedef Algebra algebra_type; typedef Algebra algebra_type;
typedef Operations operations_type; typedef Operations operations_type;
typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type; typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
typedef state_wrapper<state_type> wrapped_state_type; typedef state_wrapper<state_type> wrapped_state_type;
typedef state_wrapper<deriv_type> wrapped_deriv_type; typedef state_wrapper<deriv_type> wrapped_deriv_type;
typedef detail::adaptive_adams_coefficients<order_value, deriv_type, time_type, algebra_type, operations_type> coeff_type; typedef detail::adaptive_adams_coefficients<order_value, deriv_type, time_type, algebra_type, operations_type> coeff_type;
typedef detail::rotating_buffer<state_type, steps> error_storage_type; typedef detail::rotating_buffer<state_type, steps> error_storage_type;
typedef adaptive_adams_bashforth<order_value, state_type, value_type, state_type, value_type, algebra_type, operations_type> aab_type; typedef adaptive_adams_bashforth<order_value, state_type, value_type, state_type, value_type, algebra_type, operations_type> aab_type;
typedef adaptive_adams_moulton<order_value, state_type, value_type, state_type, value_type, algebra_type, operations_type> aam_type; typedef adaptive_adams_moulton<order_value, state_type, value_type, state_type, value_type, algebra_type, operations_type> aam_type;
typedef adaptive_adams_bashforth_moulton< Steps , State , Value , Deriv , Time, Algebra, Operations, Resizer > stepper_type; typedef adaptive_adams_bashforth_moulton< Steps , State , Value , Deriv , Time, Algebra, Operations, Resizer > stepper_type;
typedef error_stepper_tag stepper_category; typedef error_stepper_tag stepper_category;
adaptive_adams_bashforth_moulton(const algebra_type &algebra = algebra_type()) adaptive_adams_bashforth_moulton(const algebra_type &algebra = algebra_type())
:algebra_stepper_base_type( algebra ), m_aab(algebra), m_aam(algebra), :algebra_stepper_base_type( algebra ), m_aab(algebra), m_aam(algebra),
m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(), m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(),
m_coeff() m_coeff()
{}; {};
order_type order() order_type order()
{ {
return order_value + 1; return order_value + 1;
} }
order_type stepper_order() order_type stepper_order()
{ {
return order_value + 1; return order_value + 1;
} }
order_type error_order() order_type error_order()
{ {
return order_value; return order_value;
} }
template<class System> template<class System>
void do_step(System system, state_type & inOut, time_type &t, time_type &dt) void do_step(System system, state_type & inOut, time_type &t, time_type &dt)
{ {
m_xerr_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); m_xerr_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
do_step(system, inOut, t, dt, m_xerr.m_v); do_step(system, inOut, t, dt, m_xerr.m_v);
}; };
template<class System> template<class System>
void do_step(System system, const state_type & in, time_type &t, state_type & out, time_type &dt) void do_step(System system, const state_type & in, time_type &t, state_type & out, time_type &dt)
{ {
m_xerr_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); m_xerr_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
do_step(system, in, t, out, dt, m_xerr.m_v); do_step(system, in, t, out, dt, m_xerr.m_v);
}; };
template<class System> template<class System>
void do_step(System system, state_type & inOut, time_type &t, time_type &dt, state_type &xerr) 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 ) ); m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
do_step(system, inOut, t, m_xnew.m_v, dt, xerr); do_step(system, inOut, t, m_xnew.m_v, dt, xerr);
boost::numeric::odeint::copy( m_xnew.m_v , inOut); boost::numeric::odeint::copy( m_xnew.m_v , inOut);
}; };
template<class System> template<class System>
void do_step(System system, const state_type & in, time_type &t, state_type & out, time_type &dt, state_type &xerr) void do_step(System system, const state_type & in, time_type &t, state_type & out, time_type &dt, state_type &xerr)
{ {
do_step_impl(system, m_coeff, in, t, out, dt, xerr); do_step_impl(system, m_coeff, in, t, out, dt, xerr);
system(out, m_dxdt.m_v, t+dt); system(out, m_dxdt.m_v, t+dt);
m_coeff.step(m_dxdt.m_v, t+dt); m_coeff.step(m_dxdt.m_v, t+dt);
m_coeff.confirm(); m_coeff.confirm();
}; };
template<class System> template<class System>
void do_step_impl(System system, coeff_type coeff, const state_type & in, time_type &t, state_type & out, time_type &dt, state_type &xerr) void do_step_impl(System system, coeff_type coeff, const state_type & in, time_type &t, state_type & out, time_type &dt, state_type &xerr)
{ {
m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_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 ) );
if(coeff.m_effective_order == 1) if(coeff.m_effective_order == 1)
{ {
system(in, m_dxdt.m_v, t); system(in, m_dxdt.m_v, t);
coeff.step(m_dxdt.m_v, t); coeff.step(m_dxdt.m_v, t);
coeff.confirm(); coeff.confirm();
} }
// predict // predict
m_aab.do_step_impl(coeff, in, t, out, dt); m_aab.do_step_impl(coeff, in, t, out, dt);
// evaluate // evaluate
system(out, m_dxdt.m_v, t + dt); system(out, m_dxdt.m_v, t + dt);
coeff.step(m_dxdt.m_v, t + dt); coeff.step(m_dxdt.m_v, t + dt);
boost::numeric::odeint::copy( out, xerr); boost::numeric::odeint::copy( out, xerr);
// correct // correct
m_aam.do_step(coeff, in, t, out, dt); m_aam.do_step(coeff, in, t, out, dt);
// Error calculation // Error calculation
m_aab.algebra().for_each3(xerr, xerr, out, typename Operations::template scale_sum2<double, double>(-1.0, 1.0)); m_aab.algebra().for_each3(xerr, xerr, out, typename Operations::template scale_sum2<double, double>(-1.0, 1.0));
}; };
const coeff_type& coeff() const const coeff_type& coeff() const
{ {
return m_coeff; return m_coeff;
}; };
coeff_type& coeff() coeff_type& coeff()
{ {
return m_coeff; return m_coeff;
}; };
wrapped_deriv_type m_dxdt; wrapped_deriv_type m_dxdt;
private: private:
template< class StateType > template< class StateType >
bool resize_dxdt_impl( const StateType &x ) bool resize_dxdt_impl( const StateType &x )
{ {
return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() ); return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() );
}; };
template< class StateType > template< class StateType >
bool resize_xerr_impl( const StateType &x ) bool resize_xerr_impl( const StateType &x )
{ {
return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable<deriv_type>::type() ); return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable<deriv_type>::type() );
}; };
template< class StateType > template< class StateType >
bool resize_xnew_impl( const StateType &x ) bool resize_xnew_impl( const StateType &x )
{ {
return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() ); return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() );
}; };
aab_type m_aab; aab_type m_aab;
aam_type m_aam; aam_type m_aam;
resizer_type m_dxdt_resizer; resizer_type m_dxdt_resizer;
resizer_type m_xerr_resizer; resizer_type m_xerr_resizer;
resizer_type m_xnew_resizer; resizer_type m_xnew_resizer;
wrapped_state_type m_xnew; wrapped_state_type m_xnew;
wrapped_state_type m_xerr; wrapped_state_type m_xerr;
coeff_type m_coeff; coeff_type m_coeff;
}; };
}}} }}}

64
include/boost/numeric/odeint/stepper/adaptive_adams_moulton.hpp
View File

@ -28,46 +28,46 @@ class Resizer = initially_resizer
> >
class adaptive_adams_moulton: public algebra_stepper_base< Algebra , Operations > class adaptive_adams_moulton: public algebra_stepper_base< Algebra , Operations >
{ {
public: public:
static const size_t steps = Steps; static const size_t steps = Steps;
typedef unsigned short order_type; typedef unsigned short order_type;
static const order_type order_value = steps + 1; static const order_type order_value = steps + 1;
typedef State state_type; typedef State state_type;
typedef Value value_type; typedef Value value_type;
typedef Deriv deriv_type; typedef Deriv deriv_type;
typedef Time time_type; typedef Time time_type;
typedef Algebra algebra_type; typedef Algebra algebra_type;
typedef Operations operations_type; typedef Operations operations_type;
typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type; typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
typedef detail::adaptive_adams_coefficients<Steps, deriv_type, time_type, algebra_type, operations_type> coeff_type; typedef detail::adaptive_adams_coefficients<Steps, deriv_type, time_type, algebra_type, operations_type> coeff_type;
typedef adaptive_adams_moulton< Steps , State , Value , Deriv , Time , Resizer > stepper_type; typedef adaptive_adams_moulton< Steps , State , Value , Deriv , Time , Resizer > stepper_type;
adaptive_adams_moulton( const algebra_type &algebra = algebra_type()) adaptive_adams_moulton( const algebra_type &algebra = algebra_type())
:algebra_stepper_base_type( algebra ) :algebra_stepper_base_type( algebra )
{}; {};
order_type order() const order_type order() const
{ {
return order_value; return order_value;
}; };
void do_step(coeff_type & coeff, const state_type & in, void do_step(coeff_type & coeff, const state_type & in,
time_type t, state_type & out, time_type &dt) time_type t, state_type & out, time_type &dt)
{ {
coeff.poly.add_root(0); coeff.poly.add_root(0);
out = in; out = in;
// integrating // integrating
for(size_t i=0; i<coeff.m_effective_order; ++i) for(size_t i=0; i<coeff.m_effective_order; ++i)
{ {
time_type c = ((i!=coeff.m_effective_order-1)?coeff.m_c[i]:coeff.poly.evaluate_integrated(dt)); time_type c = ((i!=coeff.m_effective_order-1)?coeff.m_c[i]:coeff.poly.evaluate_integrated(dt));
this->m_algebra.for_each3(out, out, coeff.m_tss[i][coeff.m_effective_order-i-1].m_v, typename Operations::template scale_sum2<double, double>(1.0, c)); this->m_algebra.for_each3(out, out, coeff.m_tss[i][coeff.m_effective_order-i-1].m_v, typename Operations::template scale_sum2<double, double>(1.0, c));
} }
}; };
}; };
} }

154
include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp
View File

@ -18,106 +18,106 @@ namespace odeint {
template< template<
class ErrorStepper, class ErrorStepper,
class StepAdjuster = detail::pid_step_adjuster<ErrorStepper::order_value, class StepAdjuster = detail::pid_step_adjuster<ErrorStepper::order_value,
typename ErrorStepper::state_type, typename ErrorStepper::state_type,
typename ErrorStepper::time_type, typename ErrorStepper::time_type,
detail::H211PI>, detail::H211PI>,
class Resizer = initially_resizer class Resizer = initially_resizer
> >
class controlled_adams_bashforth_moulton class controlled_adams_bashforth_moulton
{ {
public: public:
typedef ErrorStepper stepper_type; typedef ErrorStepper stepper_type;
typedef typename stepper_type::state_type state_type; typedef typename stepper_type::state_type state_type;
typedef typename stepper_type::value_type value_type; typedef typename stepper_type::value_type value_type;
typedef typename stepper_type::deriv_type deriv_type; typedef typename stepper_type::deriv_type deriv_type;
typedef typename stepper_type::time_type time_type; typedef typename stepper_type::time_type time_type;
typedef typename stepper_type::algebra_type algebra_type; typedef typename stepper_type::algebra_type algebra_type;
typedef typename stepper_type::operations_type operations_type; typedef typename stepper_type::operations_type operations_type;
typedef Resizer resizer_type; typedef Resizer resizer_type;
typedef StepAdjuster step_adjuster_type; typedef StepAdjuster step_adjuster_type;
typedef controlled_stepper_tag stepper_category; typedef controlled_stepper_tag stepper_category;
typedef typename stepper_type::wrapped_state_type wrapped_state_type; typedef typename stepper_type::wrapped_state_type wrapped_state_type;
typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type; typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
typedef controlled_adams_bashforth_moulton<ErrorStepper, StepAdjuster, Resizer> controlled_stepper_type; typedef controlled_adams_bashforth_moulton<ErrorStepper, StepAdjuster, Resizer> controlled_stepper_type;
controlled_adams_bashforth_moulton(step_adjuster_type step_adjuster = step_adjuster_type()) controlled_adams_bashforth_moulton(step_adjuster_type step_adjuster = step_adjuster_type())
:m_stepper(), m_coeff(m_stepper.coeff()), :m_stepper(), m_coeff(m_stepper.coeff()),
m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(), m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(),
m_step_adjuster(step_adjuster) m_step_adjuster(step_adjuster)
{}; {};
template<class System> template<class System>
controlled_step_result try_step(System system, state_type & inOut, time_type &t, time_type &dt) 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 ) ); 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); controlled_step_result res = try_step(system, inOut, t, m_xnew.m_v, dt);
if(res == success) if(res == success)
boost::numeric::odeint::copy( m_xnew.m_v , inOut); boost::numeric::odeint::copy( m_xnew.m_v , inOut);
return res; return res;
}; };
template<class System> template<class System>
controlled_step_result try_step(System system, const state_type & in, time_type &t, state_type & out, time_type &dt) 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_xerr_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_stepper.do_step_impl(system, m_coeff, in, t, out, dt, m_xerr.m_v); m_stepper.do_step_impl(system, m_coeff, in, t, out, dt, m_xerr.m_v);
time_type dtPrev = dt; time_type dtPrev = dt;
dt = m_step_adjuster.adjust_stepsize(dt, m_xerr.m_v); dt = m_step_adjuster.adjust_stepsize(dt, m_xerr.m_v);
if(dt / dtPrev >= 0.9) if(dt / dtPrev >= 0.9)
{ {
m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_dxdt_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 ) );
system(out, m_dxdt.m_v, t+dtPrev); system(out, m_dxdt.m_v, t+dtPrev);
m_coeff.step(m_dxdt.m_v, t+dtPrev); m_coeff.step(m_dxdt.m_v, t+dtPrev);
m_coeff.confirm(); m_coeff.confirm();
t += dtPrev; t += dtPrev;
return success; return success;
} }
else else
{ {
return fail; return fail;
} }
}; };
private: private:
template< class StateType > template< class StateType >
bool resize_dxdt_impl( const StateType &x ) bool resize_dxdt_impl( const StateType &x )
{ {
return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() ); return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() );
}; };
template< class StateType > template< class StateType >
bool resize_xerr_impl( const StateType &x ) bool resize_xerr_impl( const StateType &x )
{ {
return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable<state_type>::type() ); return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable<state_type>::type() );
}; };
template< class StateType > template< class StateType >
bool resize_xnew_impl( const StateType &x ) bool resize_xnew_impl( const StateType &x )
{ {
return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() ); return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() );
}; };
stepper_type m_stepper; stepper_type m_stepper;
typename stepper_type::coeff_type &m_coeff; typename stepper_type::coeff_type &m_coeff;
wrapped_deriv_type m_dxdt; wrapped_deriv_type m_dxdt;
wrapped_state_type m_xerr; wrapped_state_type m_xerr;
wrapped_state_type m_xnew; wrapped_state_type m_xnew;
resizer_type m_dxdt_resizer; resizer_type m_dxdt_resizer;
resizer_type m_xerr_resizer; resizer_type m_xerr_resizer;
resizer_type m_xnew_resizer; resizer_type m_xnew_resizer;
step_adjuster_type m_step_adjuster; step_adjuster_type m_step_adjuster;
}; };
} }

187
include/boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp
View File

@ -32,130 +32,105 @@ class Resizer = initially_resizer
> >
struct adaptive_adams_coefficients struct adaptive_adams_coefficients
{ {
public: public:
static const size_t steps = Steps; static const size_t steps = Steps;
typedef Deriv deriv_type; typedef Deriv deriv_type;
typedef Time time_type; typedef Time time_type;
typedef state_wrapper<deriv_type> wrapped_deriv_type; typedef state_wrapper<deriv_type> wrapped_deriv_type;
typedef detail::rotating_buffer<wrapped_deriv_type, steps+1> step_storage_type; // +1 for moulton typedef detail::rotating_buffer<wrapped_deriv_type, steps+1> step_storage_type; // +1 for moulton
typedef detail::rotating_buffer<time_type, steps+1> time_storage_type; typedef detail::rotating_buffer<time_type, steps+1> time_storage_type;
typedef Algebra algebra_type; typedef Algebra algebra_type;
typedef Operations operations_type; typedef Operations operations_type;
typedef Resizer resizer_type; typedef Resizer resizer_type;
typedef adaptive_adams_coefficients<Steps, Deriv, Time, Algebra, Operations, Resizer> aac_type; typedef adaptive_adams_coefficients<Steps, Deriv, Time, Algebra, Operations, Resizer> aac_type;
typedef detail::Polynomial<steps+2, time_type> poly_type; typedef detail::Polynomial<steps+2, time_type> poly_type;
adaptive_adams_coefficients(const algebra_type &algebra = algebra_type() ) adaptive_adams_coefficients(const algebra_type &algebra = algebra_type() )
:poly(), m_effective_order(1), m_resizer(), m_algebra(algebra) :poly(), m_effective_order(1), m_resizer(), m_algebra(algebra)
{}; {};
/*aac_type &operator=(const aac_type & rhs) void step(const deriv_type &deriv, const time_type &t)
{ {
this->poly = rhs.poly; m_resizer.adjust_size( deriv , detail::bind( &aac_type::template resize_tss_impl< deriv_type > , detail::ref( *this ) , detail::_1 ) );
this->m_c = rhs.m_c;
// this->m_ss = rhs.m_ss;
for(size_t i=0; i<steps+1; ++i) m_tts[0] = t;
{ m_tss[0][0].m_v = deriv;
for(size_t j=0; j<steps+1; ++j)
{
this->m_ss[i][j] = rhs.m_ss[i][j];
}
}
this->m_ts = rhs.m_ts; for(size_t i=1; i<m_effective_order; ++i)
this->m_tss = rhs.m_tss; {
this->m_tts = rhs.m_tts; time_type dt = t - m_ts[i-1];
this->m_effective_order = rhs.m_effective_order; m_algebra.for_each3(m_tss[i][0].m_v, m_tss[i-1][0].m_v, m_ss[i-1][0].m_v, typename Operations::template scale_sum2<double, double>(1/dt, -1/dt));
}
};
void confirm()
{
for(size_t i=0; i<steps; ++i)
{
m_ss[i] = m_tss[i];
m_tss[i].rotate();
}
this->m_algebra = rhs.m_algebra; m_ts = m_tts;
this->m_resizer = rhs.m_resizer; m_tts.rotate();
return *this; if(m_effective_order < steps+1)
}*/ {
++m_effective_order;
}
};
void step(const deriv_type &deriv, const time_type &t) void reset()
{ {
m_resizer.adjust_size( deriv , detail::bind( &aac_type::template resize_tss_impl< deriv_type > , detail::ref( *this ) , detail::_1 ) ); poly.reset();
m_effective_order = 1;
};
m_tts[0] = t; void pretty_print()
m_tss[0][0].m_v = deriv; {
for(size_t k=0; k<2; ++k)
{
for(size_t i=0; i<m_effective_order; ++i)
{
for(size_t j=0; j<m_effective_order - i-1; ++j)
std::cout << m_ss[j][i].m_v[k] << " ";
std::cout << std::endl;
}
std::cout << std::endl;
}
};
for(size_t i=1; i<m_effective_order; ++i) poly_type poly;
{ time_storage_type m_c;
time_type dt = t - m_ts[i-1];
m_algebra.for_each3(m_tss[i][0].m_v, m_tss[i-1][0].m_v, m_ss[i-1][0].m_v, typename Operations::template scale_sum2<double, double>(1/dt, -1/dt));
}
};
void confirm()
{
for(size_t i=0; i<steps; ++i)
{
m_ss[i] = m_tss[i];
m_tss[i].rotate();
}
m_ts = m_tts; boost::array<step_storage_type, steps+1> m_ss;
m_tts.rotate(); time_storage_type m_ts;
if(m_effective_order < steps+1) boost::array<step_storage_type, steps+1> m_tss;
{ time_storage_type m_tts;
++m_effective_order;
}
};
void reset() size_t m_effective_order;
{
poly.reset();
m_effective_order = 1;
};
void pretty_print() private:
{ template< class StateType >
for(size_t k=0; k<2; ++k) bool resize_tss_impl( const StateType &x )
{ {
for(size_t i=0; i<m_effective_order; ++i) bool resized( false );
{ for( size_t i=0 ; i<steps+1 ; ++i )
for(size_t j=0; j<m_effective_order - i-1; ++j) {
std::cout << m_ss[j][i].m_v[k] << " "; for(size_t j=0; j<steps+1; ++j)
std::cout << std::endl; resized |= adjust_size_by_resizeability( m_tss[i][j] , x , typename is_resizeable<deriv_type>::type() );
}
std::cout << std::endl; m_tts[i] = 0;
} }
}; return resized;
};
poly_type poly; resizer_type m_resizer;
time_storage_type m_c; algebra_type m_algebra;
boost::array<step_storage_type, steps+1> m_ss;
time_storage_type m_ts;
boost::array<step_storage_type, steps+1> m_tss;
time_storage_type m_tts;
size_t m_effective_order;
private:
template< class StateType >
bool resize_tss_impl( const StateType &x )
{
bool resized( false );
for( size_t i=0 ; i<steps+1 ; ++i )
{
for(size_t j=0; j<steps+1; ++j)
resized |= adjust_size_by_resizeability( m_tss[i][j] , x , typename is_resizeable<deriv_type>::type() );
m_tts[i] = 0;
}
return resized;
};
resizer_type m_resizer;
algebra_type m_algebra;
}; };
} }

116
include/boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp
View File

@ -19,79 +19,79 @@ size_t Type = H211PI
> >
struct pid_step_adjuster struct pid_step_adjuster
{ {
public: public:
const static size_t steps = Steps; const static size_t steps = Steps;
typedef State state_type; typedef State state_type;
typedef Time time_type; typedef Time time_type;
typedef rotating_buffer<state_type, 3> error_storage_type; typedef rotating_buffer<state_type, 3> error_storage_type;
typedef rotating_buffer<time_type, 3> time_storage_type; typedef rotating_buffer<time_type, 3> time_storage_type;
typedef pid_step_adjuster_coefficients<Type> coeff_type; typedef pid_step_adjuster_coefficients<Type> coeff_type;
pid_step_adjuster(double tol = 1e-5, time_type dtmax = 1.0) pid_step_adjuster(double tol = 1e-5, time_type dtmax = 1.0)
:m_dtmax(dtmax), m_error_storage(), m_time_storage(), m_init(0), m_tol(tol) :m_dtmax(dtmax), m_error_storage(), m_time_storage(), m_init(0), m_tol(tol)
{}; {};
time_type adjust_stepsize(time_type dt, const state_type &err) time_type adjust_stepsize(time_type dt, const state_type &err)
{ {
m_error_storage[0] = err; m_error_storage[0] = err;
m_time_storage[0] = dt; m_time_storage[0] = dt;
double ratio = 100; double ratio = 100;
double r; double r;
for(size_t i=0; i<m_error_storage[0].size(); ++i) for(size_t i=0; i<m_error_storage[0].size(); ++i)
{ {
if(m_init >= 2) if(m_init >= 2)
{ {
r = pow(fabs(m_tol/m_error_storage[0][i]), m_coeff[0]/(steps + 1)) * r = pow(fabs(m_tol/m_error_storage[0][i]), m_coeff[0]/(steps + 1)) *
pow(fabs(m_tol/m_error_storage[1][i]), m_coeff[1]/(steps + 1)) * pow(fabs(m_tol/m_error_storage[1][i]), m_coeff[1]/(steps + 1)) *
pow(fabs(m_tol/m_error_storage[2][i]), m_coeff[2]/(steps + 1)) * pow(fabs(m_tol/m_error_storage[2][i]), m_coeff[2]/(steps + 1)) *
pow(fabs(m_time_storage[0]/m_time_storage[1]), -m_coeff[3]/(steps + 1))* pow(fabs(m_time_storage[0]/m_time_storage[1]), -m_coeff[3]/(steps + 1))*
pow(fabs(m_time_storage[1]/m_time_storage[2]), -m_coeff[4]/(steps + 1)); pow(fabs(m_time_storage[1]/m_time_storage[2]), -m_coeff[4]/(steps + 1));
} }
else else
{ {
r = pow(fabs(m_tol/m_error_storage[0][i]), 0.7/(steps + 1)); // purely integrating controller for startup r = pow(fabs(m_tol/m_error_storage[0][i]), 0.7/(steps + 1)); // purely integrating controller for startup
} }
if(r<ratio) if(r<ratio)
ratio = r; ratio = r;
} }
double kappa = 1.0; double kappa = 1.0;
ratio = 1.0 + kappa*atan((ratio - 1)/kappa); ratio = 1.0 + kappa*atan((ratio - 1)/kappa);
if(ratio*dt >= m_dtmax) if(ratio*dt >= m_dtmax)
{ {
dt = m_dtmax; dt = m_dtmax;
} }
if(ratio >= 0.9) if(ratio >= 0.9)
{ {
m_error_storage.rotate(); m_error_storage.rotate();
m_time_storage.rotate(); m_time_storage.rotate();
++m_init; ++m_init;
} }
else else
{ {
m_init = 0; m_init = 0;
} }
return dt*ratio; return dt*ratio;
}; };
private: private:
time_type m_dtmax; time_type m_dtmax;
error_storage_type m_error_storage; error_storage_type m_error_storage;
time_storage_type m_time_storage; time_storage_type m_time_storage;
size_t m_init; size_t m_init;
double m_tol; double m_tol;
coeff_type m_coeff; coeff_type m_coeff;
}; };
} }

180
include/boost/numeric/odeint/stepper/detail/pid_step_adjuster_coefficients.hpp
View File

@ -9,15 +9,15 @@ namespace odeint{
namespace detail{ namespace detail{
enum adjuster_type{ enum adjuster_type{
BASIC, BASIC,
H0211, H0211,
H211b, H211b,
H211PI, H211PI,
H0312, H0312,
H312b, H312b,
H312PID, H312PID,
H0321, H0321,
H321 H321
}; };
template<int Type> template<int Type>
@ -27,135 +27,135 @@ template<>
class pid_step_adjuster_coefficients<BASIC> : public boost::array<double, 5> class pid_step_adjuster_coefficients<BASIC> : public boost::array<double, 5>
{ {
public: public:
pid_step_adjuster_coefficients() pid_step_adjuster_coefficients()
: boost::array<double, 5>() : boost::array<double, 5>()
{ {
(*this)[0] = 1; (*this)[0] = 1;
(*this)[1] = 0; (*this)[1] = 0;
(*this)[2] = 0; (*this)[2] = 0;
(*this)[3] = 0; (*this)[3] = 0;
(*this)[4] = 0; (*this)[4] = 0;
} }
}; };
template<> template<>
class pid_step_adjuster_coefficients<H0211> : public boost::array<double, 5> class pid_step_adjuster_coefficients<H0211> : public boost::array<double, 5>
{ {
public: public:
pid_step_adjuster_coefficients() pid_step_adjuster_coefficients()
: boost::array<double, 5>() : boost::array<double, 5>()
{ {
(*this)[0] = 1.0/2.0; (*this)[0] = 1.0/2.0;
(*this)[1] = 1.0/2.0; (*this)[1] = 1.0/2.0;
(*this)[2] = 0; (*this)[2] = 0;
(*this)[3] = 1.0/2.0; (*this)[3] = 1.0/2.0;
(*this)[4] = 0; (*this)[4] = 0;
} }
}; };
template<> template<>
class pid_step_adjuster_coefficients<H211b> : public boost::array<double, 5> class pid_step_adjuster_coefficients<H211b> : public boost::array<double, 5>
{ {
public: public:
pid_step_adjuster_coefficients() pid_step_adjuster_coefficients()
: boost::array<double, 5>() : boost::array<double, 5>()
{ {
(*this)[0] = 1.0/5.0; (*this)[0] = 1.0/5.0;
(*this)[1] = 2.0/5.0; (*this)[1] = 2.0/5.0;
(*this)[2] = 0; (*this)[2] = 0;
(*this)[3] = 1.0/5.0; (*this)[3] = 1.0/5.0;
(*this)[4] = 0; (*this)[4] = 0;
} }
}; };
template<> template<>
class pid_step_adjuster_coefficients<H211PI> : public boost::array<double, 5> class pid_step_adjuster_coefficients<H211PI> : public boost::array<double, 5>
{ {
public: public:
pid_step_adjuster_coefficients() pid_step_adjuster_coefficients()
: boost::array<double, 5>() : boost::array<double, 5>()
{ {
(*this)[0] = 1.0/6.0; (*this)[0] = 1.0/6.0;
(*this)[1] = 2.0/6.0; (*this)[1] = 2.0/6.0;
(*this)[2] = 0; (*this)[2] = 0;
(*this)[3] = 0; (*this)[3] = 0;
(*this)[4] = 0; (*this)[4] = 0;
} }
}; };
template<> template<>
class pid_step_adjuster_coefficients<H0312> : public boost::array<double, 5> class pid_step_adjuster_coefficients<H0312> : public boost::array<double, 5>
{ {
public: public:
pid_step_adjuster_coefficients() pid_step_adjuster_coefficients()
: boost::array<double, 5>() : boost::array<double, 5>()
{ {
(*this)[0] = 1.0/4.0; (*this)[0] = 1.0/4.0;
(*this)[1] = 2.0/2.0; (*this)[1] = 2.0/2.0;
(*this)[2] = 1.0/4.0; (*this)[2] = 1.0/4.0;
(*this)[3] = 3.0/4.0; (*this)[3] = 3.0/4.0;
(*this)[4] = 1.0/4.0; (*this)[4] = 1.0/4.0;
} }
}; };
template<> template<>
class pid_step_adjuster_coefficients<H312b> : public boost::array<double, 5> class pid_step_adjuster_coefficients<H312b> : public boost::array<double, 5>
{ {
public: public:
pid_step_adjuster_coefficients() pid_step_adjuster_coefficients()
: boost::array<double, 5>() : boost::array<double, 5>()
{ {
(*this)[0] = 1.0/6.0; (*this)[0] = 1.0/6.0;
(*this)[1] = 2.0/6.0; (*this)[1] = 2.0/6.0;
(*this)[2] = 1.0/6.0; (*this)[2] = 1.0/6.0;
(*this)[3] = 3.0/6.0; (*this)[3] = 3.0/6.0;
(*this)[4] = 1.0/6.0; (*this)[4] = 1.0/6.0;
} }
}; };
template<> template<>
class pid_step_adjuster_coefficients<H312PID> : public boost::array<double, 5> class pid_step_adjuster_coefficients<H312PID> : public boost::array<double, 5>
{ {
public: public:
pid_step_adjuster_coefficients() pid_step_adjuster_coefficients()
: boost::array<double, 5>() : boost::array<double, 5>()
{ {
(*this)[0] = 1.0/18.0; (*this)[0] = 1.0/18.0;
(*this)[1] = 2.0/9.0; (*this)[1] = 2.0/9.0;
(*this)[2] = 1.0/18.0; (*this)[2] = 1.0/18.0;
(*this)[3] = 0; (*this)[3] = 0;
(*this)[4] = 0; (*this)[4] = 0;
} }
}; };
template<> template<>
class pid_step_adjuster_coefficients<H0321> : public boost::array<double, 5> class pid_step_adjuster_coefficients<H0321> : public boost::array<double, 5>
{ {
public: public:
pid_step_adjuster_coefficients() pid_step_adjuster_coefficients()
: boost::array<double, 5>() : boost::array<double, 5>()
{ {
(*this)[0] = 5.0/4.0; (*this)[0] = 5.0/4.0;
(*this)[1] = 1.0/2.0; (*this)[1] = 1.0/2.0;
(*this)[2] = -3.0/4.0; (*this)[2] = -3.0/4.0;
(*this)[3] = -1.0/4.0; (*this)[3] = -1.0/4.0;
(*this)[4] = -3.0/4.0; (*this)[4] = -3.0/4.0;
} }
}; };
template<> template<>
class pid_step_adjuster_coefficients<H321> : public boost::array<double, 5> class pid_step_adjuster_coefficients<H321> : public boost::array<double, 5>
{ {
public: public:
pid_step_adjuster_coefficients() pid_step_adjuster_coefficients()
: boost::array<double, 5>() : boost::array<double, 5>()
{ {
(*this)[0] = 1.0/3.0; (*this)[0] = 1.0/3.0;
(*this)[1] = 1.0/18.0; (*this)[1] = 1.0/18.0;
(*this)[2] = -5.0/18.0; (*this)[2] = -5.0/18.0;
(*this)[3] = -5.0/16.0; (*this)[3] = -5.0/16.0;
(*this)[4] = -1.0/6.0; (*this)[4] = -1.0/6.0;
} }
}; };
} }

206
include/boost/numeric/odeint/stepper/detail/polynomial.hpp
View File

@ -15,127 +15,127 @@ class Time
> >
class Polynomial class Polynomial
{ {
public: public:
typedef Time time_type; typedef Time time_type;
typedef Polynomial<order, Time> poly_type; typedef Polynomial<order, Time> poly_type;
Polynomial() Polynomial()
{ {
for(size_t i = 0; i<order; ++i) for(size_t i = 0; i<order; ++i)
m_coeff[i] = 0; m_coeff[i] = 0;
// start with 'constant' polynomial // start with 'constant' polynomial
m_coeff[order-1] = 1; m_coeff[order-1] = 1;
}; };
Polynomial(boost::array<time_type, order> coeff) Polynomial(boost::array<time_type, order> coeff)
{ {
m_coeff = coeff; m_coeff = coeff;
}; };
Polynomial<order+1, time_type> integrate() Polynomial<order+1, time_type> integrate()
{ {
boost::array<time_type, order + 1> coeff_int; boost::array<time_type, order + 1> coeff_int;
coeff_int[order] = 0; coeff_int[order] = 0;
for(size_t i=0; i<order; ++i) for(size_t i=0; i<order; ++i)
{ {
if(i != order) if(i != order)
coeff_int[i] = m_coeff[i] / (order - i); coeff_int[i] = m_coeff[i] / (order - i);
} }
Polynomial<order+1, time_type> polyInt(coeff_int); Polynomial<order+1, time_type> polyInt(coeff_int);
return polyInt; return polyInt;
}; };
time_type evaluate(const time_type &t) time_type evaluate(const time_type &t)
{ {
// fma: x*y+z // fma: x*y+z
m_res = m_coeff[0]; m_res = m_coeff[0];
for(size_t i=1; i<order; ++i) for(size_t i=1; i<order; ++i)
{ {
// m_res = fma(m_res, t, m_coeff[i]); // m_res = fma(m_res, t, m_coeff[i]);
m_res = m_coeff[i] + m_res * t; m_res = m_coeff[i] + m_res * t;
} }
return m_res; return m_res;
}; };
time_type evaluate_integrated(const time_type &t) time_type evaluate_integrated(const time_type &t)
{ {
m_res = m_coeff[0]/order; m_res = m_coeff[0]/order;
for(size_t i=1; i<order+1; ++i) for(size_t i=1; i<order+1; ++i)
{ {
// m_res = fma(m_res, t, ((i >= order)?0:m_coeff[i]/(order-i))); // m_res = fma(m_res, t, ((i >= order)?0:m_coeff[i]/(order-i)));
m_res = ((i >= order)?0:m_coeff[i]/(order-i)) + m_res * t; m_res = ((i >= order)?0:m_coeff[i]/(order-i)) + m_res * t;
} }
return m_res; return m_res;
} }
void add_root(const time_type &root) void add_root(const time_type &root)
{ {
for(size_t j=0; j<order; ++j) // updating all coefficients for(size_t j=0; j<order; ++j) // updating all coefficients
{ {
m_coeff[j] = -m_coeff[j]*root + ((j<order-1)?m_coeff[j+1]:0); m_coeff[j] = -m_coeff[j]*root + ((j<order-1)?m_coeff[j+1]:0);
} }
}; };
void reset() void reset()
{ {
for(size_t i = 0; i<order; ++i) for(size_t i = 0; i<order; ++i)
m_coeff[i] = 0; m_coeff[i] = 0;
m_coeff[order-1] = 1; m_coeff[order-1] = 1;
}; };
void remove_root(const time_type &root) void remove_root(const time_type &root)
{ {
// 'reverse' add_root; synthetic division // 'reverse' add_root; synthetic division
time_type tmp[2]; time_type tmp[2];
m_coeff[0] = 0; m_coeff[0] = 0;
for(size_t j=0; j<order; ++j) for(size_t j=0; j<order; ++j)
{ {
tmp[0] = (j!=0)?tmp[1]:0; tmp[0] = (j!=0)?tmp[1]:0;
tmp[1] = m_coeff[j+1]; tmp[1] = m_coeff[j+1];
m_coeff[j+1] = tmp[0] + root * m_coeff[j]; m_coeff[j+1] = tmp[0] + root * m_coeff[j];
} }
}; };
void pretty_print() void pretty_print()
{ {
for(size_t i=0; i<order; ++i) for(size_t i=0; i<order; ++i)
{ {
std::cout << m_coeff[i]; std::cout << m_coeff[i];
std::cout << "x^" << order - i - 1; std::cout << "x^" << order - i - 1;
if(i != order - 1) if(i != order - 1)
std::cout<< " + "; std::cout<< " + ";
else else
std::cout << std::endl; std::cout << std::endl;
} }
}; };
static Polynomial<order, time_type> from_roots(time_type * roots, unsigned short num_roots) static Polynomial<order, time_type> from_roots(time_type * roots, unsigned short num_roots)
{ {
time_type coeff[order] = {0}; time_type coeff[order] = {0};
coeff[order-1] = 1; coeff[order-1] = 1;
for(size_t i=0; i<num_roots; ++i) // going over all roots for(size_t i=0; i<num_roots; ++i) // going over all roots
{ {
for(size_t j=0; j<order; ++j) // updating all coefficients for(size_t j=0; j<order; ++j) // updating all coefficients
{ {
coeff[j] = -coeff[j]*roots[i] + ((j<order)?coeff[j+1]:0); coeff[j] = -coeff[j]*roots[i] + ((j<order)?coeff[j+1]:0);
} }
} }
Polynomial<order, time_type> poly(coeff); Polynomial<order, time_type> poly(coeff);
return poly; return poly;
}; };
// first element is highest order // first element is highest order
boost::array<time_type, order> m_coeff; boost::array<time_type, order> m_coeff;
private: private:
time_type m_res; time_type m_res;
}; };
}}}} }}}}

48
test/adaptive_adams_bashforth.cpp
View File

@ -27,33 +27,33 @@ BOOST_AUTO_TEST_SUITE( adaptive_adams_bashforth_test )
BOOST_AUTO_TEST_CASE( test_instantiation ) BOOST_AUTO_TEST_CASE( test_instantiation )
{ {
typedef double time_type; typedef double time_type;
typedef std::vector<double> state_type; typedef std::vector<double> state_type;
adaptive_adams_bashforth<1, state_type> s1; adaptive_adams_bashforth<1, state_type> s1;
adaptive_adams_bashforth<2, state_type> s2; adaptive_adams_bashforth<2, state_type> s2;
adaptive_adams_bashforth<3, state_type> s3; adaptive_adams_bashforth<3, state_type> s3;
adaptive_adams_bashforth<4, state_type> s4; adaptive_adams_bashforth<4, state_type> s4;
adaptive_adams_bashforth<5, state_type> s5; adaptive_adams_bashforth<5, state_type> s5;
adaptive_adams_bashforth<6, state_type> s6; adaptive_adams_bashforth<6, state_type> s6;
adaptive_adams_bashforth<7, state_type> s7; adaptive_adams_bashforth<7, state_type> s7;
adaptive_adams_bashforth<8, state_type> s8; adaptive_adams_bashforth<8, state_type> s8;
adaptive_adams_bashforth<9, state_type> s9; adaptive_adams_bashforth<9, state_type> s9;
state_type x0; state_type x0;
x0.push_back(0); x0.push_back(0);
time_type t0 = 0.0; time_type t0 = 0.0;
time_type dt = 0.1; time_type dt = 0.1;
s1.do_step(const_sys(), x0, t0, dt); s1.do_step(const_sys(), x0, t0, dt);
s2.do_step(const_sys(), x0, t0, dt); s2.do_step(const_sys(), x0, t0, dt);
s3.do_step(const_sys(), x0, t0, dt); s3.do_step(const_sys(), x0, t0, dt);
s4.do_step(const_sys(), x0, t0, dt); s4.do_step(const_sys(), x0, t0, dt);
s5.do_step(const_sys(), x0, t0, dt); s5.do_step(const_sys(), x0, t0, dt);
s6.do_step(const_sys(), x0, t0, dt); s6.do_step(const_sys(), x0, t0, dt);
s7.do_step(const_sys(), x0, t0, dt); s7.do_step(const_sys(), x0, t0, dt);
s8.do_step(const_sys(), x0, t0, dt); s8.do_step(const_sys(), x0, t0, dt);
s9.do_step(const_sys(), x0, t0, dt); s9.do_step(const_sys(), x0, t0, dt);
} }
BOOST_AUTO_TEST_CASE( test_copy ) BOOST_AUTO_TEST_CASE( test_copy )

70
test/adaptive_adams_coefficients.cpp
View File

@ -25,58 +25,58 @@ BOOST_AUTO_TEST_SUITE( adaptive_adams_coefficients_test )
typedef boost::mpl::range_c< size_t , 2 , 10 > vector_of_steps; typedef boost::mpl::range_c< size_t , 2 , 10 > vector_of_steps;
BOOST_AUTO_TEST_CASE_TEMPLATE( test_step, step_type, vector_of_steps ) BOOST_AUTO_TEST_CASE_TEMPLATE( test_step, step_type, vector_of_steps )
{ {
const static size_t steps = step_type::value; const static size_t steps = step_type::value;
typedef std::vector<double> deriv_type; typedef std::vector<double> deriv_type;
typedef double time_type; typedef double time_type;
typedef detail::adaptive_adams_coefficients<steps, deriv_type, time_type> aac_type; typedef detail::adaptive_adams_coefficients<steps, deriv_type, time_type> aac_type;
std::vector<double> deriv; std::vector<double> deriv;
deriv.push_back(-1); deriv.push_back(-1);
aac_type coeff; aac_type coeff;
for(size_t i=0; i<steps; ++i) for(size_t i=0; i<steps; ++i)
{ {
coeff.step(deriv, i); coeff.step(deriv, i);
BOOST_CHECK_EQUAL(coeff.m_tts[0], i); BOOST_CHECK_EQUAL(coeff.m_tts[0], i);
coeff.confirm(); coeff.confirm();
BOOST_CHECK_EQUAL(coeff.m_ts[0], i); BOOST_CHECK_EQUAL(coeff.m_ts[0], i);
} }
BOOST_CHECK_EQUAL(coeff.m_ss[0][0].m_v[0], -1); BOOST_CHECK_EQUAL(coeff.m_ss[0][0].m_v[0], -1);
BOOST_CHECK_EQUAL(coeff.m_ss[1][0].m_v[0], 0); BOOST_CHECK_EQUAL(coeff.m_ss[1][0].m_v[0], 0);
} }
BOOST_AUTO_TEST_CASE( test_copy ) BOOST_AUTO_TEST_CASE( test_copy )
{ {
typedef std::vector<double> deriv_type; typedef std::vector<double> deriv_type;
typedef double time_type; typedef double time_type;
typedef detail::adaptive_adams_coefficients<3, deriv_type, time_type> aac_type; typedef detail::adaptive_adams_coefficients<3, deriv_type, time_type> aac_type;
aac_type c1; aac_type c1;
deriv_type deriv(1); deriv_type deriv(1);
deriv[0] = 1.0; deriv[0] = 1.0;
c1.step(deriv, 0.0); c1.step(deriv, 0.0);
c1.confirm(); c1.confirm();
c1.step(deriv, 1.0); c1.step(deriv, 1.0);
c1.confirm(); c1.confirm();
c1.step(deriv, 2.0); c1.step(deriv, 2.0);
c1.confirm(); c1.confirm();
aac_type c2(c1); aac_type c2(c1);
BOOST_CHECK_EQUAL(c1.m_ss[0][0].m_v[0], c2.m_ss[0][0].m_v[0]); BOOST_CHECK_EQUAL(c1.m_ss[0][0].m_v[0], c2.m_ss[0][0].m_v[0]);
BOOST_CHECK(&(c1.m_ss[0][0].m_v) != &(c2.m_ss[0][0].m_v)); BOOST_CHECK(&(c1.m_ss[0][0].m_v) != &(c2.m_ss[0][0].m_v));
aac_type c3; aac_type c3;
deriv_type *p1 = &(c3.m_ss[0][0].m_v); deriv_type *p1 = &(c3.m_ss[0][0].m_v);
c3 = c1; c3 = c1;
// BOOST_CHECK(p1 == (&(c3.m_ss[0][0].m_v))); // BOOST_CHECK(p1 == (&(c3.m_ss[0][0].m_v)));
BOOST_CHECK_EQUAL(c1.m_ss[0][0].m_v[0], c3.m_ss[0][0].m_v[0]); BOOST_CHECK_EQUAL(c1.m_ss[0][0].m_v[0], c3.m_ss[0][0].m_v[0]);
} }
BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE_END()

38
test/controlled_adams_bashforth_moulton.cpp
View File

@ -29,28 +29,28 @@ BOOST_AUTO_TEST_SUITE( controlled_adams_bashforth_moulton_test )
BOOST_AUTO_TEST_CASE( test_instantiation ) BOOST_AUTO_TEST_CASE( test_instantiation )
{ {
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<1, state_type> > s1; controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<1, state_type> > s1;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<2, state_type> > s2; controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<2, state_type> > s2;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<3, state_type> > s3; controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<3, state_type> > s3;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<4, state_type> > s4; controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<4, state_type> > s4;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<5, state_type> > s5; controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<5, state_type> > s5;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<6, state_type> > s6; controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<6, state_type> > s6;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<7, state_type> > s7; controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<7, state_type> > s7;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<8, state_type> > s8; controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<8, state_type> > s8;
controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<9, state_type> > s9; controlled_adams_bashforth_moulton<adaptive_adams_bashforth_moulton<9, state_type> > s9;
state_type x = {{ 10.0 }}; state_type x = {{ 10.0 }};
value_type t = 0.0 , dt = 0.01; value_type t = 0.0 , dt = 0.01;
s1.try_step(const_sys(), x, t, dt); s1.try_step(const_sys(), x, t, dt);
s2.try_step(const_sys(), x, t, dt); s2.try_step(const_sys(), x, t, dt);
s3.try_step(const_sys(), x, t, dt); s3.try_step(const_sys(), x, t, dt);
s4.try_step(const_sys(), x, t, dt); s4.try_step(const_sys(), x, t, dt);
s5.try_step(const_sys(), x, t, dt); s5.try_step(const_sys(), x, t, dt);
s6.try_step(const_sys(), x, t, dt); s6.try_step(const_sys(), x, t, dt);
s7.try_step(const_sys(), x, t, dt); s7.try_step(const_sys(), x, t, dt);
s8.try_step(const_sys(), x, t, dt); s8.try_step(const_sys(), x, t, dt);
s9.try_step(const_sys(), x, t, dt); s9.try_step(const_sys(), x, t, dt);
} }
BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE_END()

62
test/polynomial.cpp
View File

@ -16,61 +16,61 @@ BOOST_AUTO_TEST_SUITE( polynomial_test )
BOOST_AUTO_TEST_CASE( test_init ) BOOST_AUTO_TEST_CASE( test_init )
{ {
Polynomial<3, double> p1; Polynomial<3, double> p1;
boost::array<double, 3> coeff = {0, 1, 2}; boost::array<double, 3> coeff = {0, 1, 2};
Polynomial<3, double> p2(coeff); Polynomial<3, double> p2(coeff);
} }
BOOST_AUTO_TEST_CASE( test_add_roots ) BOOST_AUTO_TEST_CASE( test_add_roots )
{ {
Polynomial<3, double> poly; Polynomial<3, double> poly;
poly.add_root(1); poly.add_root(1);
poly.add_root(0); poly.add_root(0);
} }
BOOST_AUTO_TEST_CASE( test_copy ) BOOST_AUTO_TEST_CASE( test_copy )
{ {
typedef Polynomial<5, double> poly_type; typedef Polynomial<5, double> poly_type;
poly_type p1; poly_type p1;
p1.add_root(1); p1.add_root(1);
p1.add_root(0); p1.add_root(0);
poly_type p2(p1); poly_type p2(p1);
BOOST_CHECK_EQUAL(p1.m_coeff[0], p2.m_coeff[0]); BOOST_CHECK_EQUAL(p1.m_coeff[0], p2.m_coeff[0]);
BOOST_CHECK_EQUAL(p1.m_coeff[1], p2.m_coeff[1]); BOOST_CHECK_EQUAL(p1.m_coeff[1], p2.m_coeff[1]);
poly_type p3; poly_type p3;
double* a1 = &(p3.m_coeff[0]); double* a1 = &(p3.m_coeff[0]);
p3 = p1; p3 = p1;
BOOST_CHECK(a1 == &(p3.m_coeff[0])); BOOST_CHECK(a1 == &(p3.m_coeff[0]));
BOOST_CHECK_EQUAL(p1.m_coeff[0], p3.m_coeff[0]); BOOST_CHECK_EQUAL(p1.m_coeff[0], p3.m_coeff[0]);
BOOST_CHECK_EQUAL(p1.m_coeff[1], p3.m_coeff[1]); BOOST_CHECK_EQUAL(p1.m_coeff[1], p3.m_coeff[1]);
} }
BOOST_AUTO_TEST_CASE( test_remove ) BOOST_AUTO_TEST_CASE( test_remove )
{ {
} }
BOOST_AUTO_TEST_CASE( test_integrate ) BOOST_AUTO_TEST_CASE( test_integrate )
{ {
} }
BOOST_AUTO_TEST_CASE( test_evaluate ) BOOST_AUTO_TEST_CASE( test_evaluate )
{ {
boost::array<double, 3> c1 = {2, 1, 0}; boost::array<double, 3> c1 = {2, 1, 0};
Polynomial<3, double> p1(c1); Polynomial<3, double> p1(c1);
BOOST_CHECK_EQUAL(p1.evaluate(0), 0); BOOST_CHECK_EQUAL(p1.evaluate(0), 0);
BOOST_CHECK_EQUAL(p1.evaluate(1), 3); BOOST_CHECK_EQUAL(p1.evaluate(1), 3);
BOOST_CHECK_EQUAL(p1.evaluate(2), 10); BOOST_CHECK_EQUAL(p1.evaluate(2), 10);
boost::array<double, 5> c2 = {0.001, 10, 0, 3, 1}; boost::array<double, 5> c2 = {0.001, 10, 0, 3, 1};
Polynomial<5, double> p2(c2); Polynomial<5, double> p2(c2);
BOOST_CHECK_CLOSE(p2.evaluate(0), 1, 0.001); BOOST_CHECK_CLOSE(p2.evaluate(0), 1, 0.001);
BOOST_CHECK_CLOSE(p2.evaluate(0.001), 1.003, 0.001); BOOST_CHECK_CLOSE(p2.evaluate(0.001), 1.003, 0.001);
BOOST_CHECK_CLOSE(p2.evaluate(1.001), 14.034, 0.001); BOOST_CHECK_CLOSE(p2.evaluate(1.001), 14.034, 0.001);
} }
BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE_END()