boost/odeint
2
0
Fork 0
mirror of https://github.com/boostorg/odeint.git synced 2026-01-19 04:22:12 +00:00
Code Issues Projects Releases Wiki Activity

moving to new formulation

Browse Source
This commit is contained in:
Valentin Hartmann
2017-07-02 15:14:09 +02:00
parent 7fe4477acf
commit c5aab45386
7 changed files with 318 additions and 602 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
include/boost/numeric/odeint.hpp
View File

@@ -51,7 +51,6 @@
#include <boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/stepper/adaptive_adams_bashforth.hpp>
#include <boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp>
#include <boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp>

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

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

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

@@ -1,44 +1,46 @@
#ifndef STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
#define STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_ADAPTIVE_ADAMS_BASHFORTH_MOULTON_HPP_INCLUDED
#include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
#include <boost/numeric/odeint/stepper/adaptive_adams_bashforth.hpp>
#include <boost/numeric/odeint/stepper/adaptive_adams_moulton.hpp>
#include <boost/numeric/odeint/stepper/detail/adaptive_adams_coefficients.hpp>
#include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
#include <boost/numeric/odeint/util/bind.hpp>
#include <boost/numeric/odeint/util/copy.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/stepper/controlled_step_result.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/algebra/algebra_dispatcher.hpp>
#include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/util/bind.hpp>
#include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
#include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
#include <iostream>
namespace boost {
namespace numeric {
namespace odeint {
template<
size_t Steps,
class State,
class Value = double,
class Deriv = State,
class Time = Value,
class Algebra = typename algebra_dispatcher< State >::algebra_type,
class Algebra = typename algebra_dispatcher< State >::algebra_type ,
class Operations = typename operations_dispatcher< State >::operations_type ,
class Resizer = initially_resizer
class Resizer = initially_resizer
>
class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , Operations >
{
public:
static const size_t steps = Steps;
typedef unsigned short order_type;
static const order_type order_value = steps;
@@ -46,132 +48,143 @@ public:
typedef Value value_type;
typedef Deriv deriv_type;
typedef Time time_type;
typedef Resizer resizer_type;
typedef Algebra algebra_type;
typedef Operations operations_type;
typedef state_wrapper< state_type > wrapped_state_type;
typedef state_wrapper< deriv_type > wrapped_deriv_type;
typedef boost::array<wrapped_state_type, 3> error_storage_type;
typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
typedef state_wrapper<state_type> wrapped_state_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::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_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 typename algebra_stepper_base_type::algebra_type algebra_type;
typedef typename algebra_stepper_base_type::operations_type operations_type;
typedef Resizer resizer_type;
typedef error_stepper_tag stepper_category;
adaptive_adams_bashforth_moulton(const algebra_type &algebra = algebra_type())
:algebra_stepper_base_type( algebra ), m_aab(algebra), m_aam(algebra),
m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(),
m_coeff()
typedef detail::adaptive_adams_coefficients< Steps, Deriv, Value, Time, Algebra, Operations, Resizer> coeff_type;
typedef adaptive_adams_bashforth_moulton< Steps , State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_type;
order_type order() const { return order_value; };
order_type stepper_order() const {return order_value + 1; };
order_type error_order() const {return order_value; };
adaptive_adams_bashforth_moulton( const algebra_type &algebra = algebra_type() )
:algebra_stepper_base_type( algebra ), m_coeff(),
m_dxdt_resizer(), m_xnew_resizer(), m_xerr_resizer()
{};
order_type order()
{
return order_value;
}
order_type stepper_order()
{
return order_value + 1;
}
order_type error_order()
{
return order_value;
}
template<class System>
void do_step(System system, state_type & inOut, time_type t, time_type &dt)
{
m_xerr_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
do_step(system, inOut, t, dt, m_xerr.m_v);
};
template<class System>
void do_step(System system, const state_type & in, time_type t, state_type & out, time_type &dt)
{
m_xerr_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
do_step(system, in, t, out, dt, m_xerr.m_v);
};
template<class System>
void do_step(System system, state_type & inOut, time_type t, time_type &dt, state_type &xerr)
template< class System >
void do_step(System system, state_type &inOut, time_type t, time_type dt )
{
m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
do_step(system, inOut, t, m_xnew.m_v, dt, xerr);
do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr[1].m_v);
boost::numeric::odeint::copy( m_xnew.m_v , inOut);
};
template<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, m_coeff, in, t, out, dt, xerr);
system(out, m_dxdt.m_v, t+dt);
m_coeff.step(m_dxdt.m_v, t+dt);
m_coeff.confirm();
template< class System >
void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt )
{
do_step(system, in, t, out, dt, m_xerr);
};
template<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)
template< class System >
void do_step(System system, state_type &inOut, time_type t, time_type dt, state_type &xerr)
{
m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr[1].m_v);
boost::numeric::odeint::copy( m_xnew.m_v , inOut);
boost::numeric::odeint::copy( m_xerr[1].m_v , xerr);
};
template< class System >
void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt , state_type &xerr)
{
do_step_impl(system, in, t, out, dt, m_xerr);
boost::numeric::odeint::copy( m_xerr[1].m_v , xerr);
if(coeff.m_effective_order == 1)
system(out, m_dxdt.m_v, t+dt);
m_coeff.do_step(m_dxdt.m_v);
m_coeff.confirm();
if(m_coeff.m_eo < order_value)
m_coeff.m_eo ++;
};
template<class ExplicitStepper, class System>
void initialize(ExplicitStepper stepper, System system, state_type &inOut, time_type &t, time_type dt)
{
m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
for( size_t i=0 ; i<order_value ; ++i )
{
system(in, m_dxdt.m_v, t);
sys( inOut , m_xnew.m_v , t );
stepper.do_step_dxdt_impl( system, inOut, m_xnew.m_v, t, dt );
m_coeff.predict(t, dt);
m_coeff.do_step(m_dxdt.m_v);
m_coeff.confirm();
coeff.step(m_dxdt.m_v, t);
coeff.confirm();
t += dt;
}
// predict
m_aab.do_step_impl(coeff, in, t, out, dt);
// evaluate
system(out, m_dxdt.m_v, t + dt);
coeff.step(m_dxdt.m_v, t + dt);
boost::numeric::odeint::copy( out, xerr);
// correct
m_aam.do_step(coeff, in, t, out, dt);
// Error calculation
m_aab.algebra().for_each3(xerr, xerr, out, typename Operations::template scale_sum2<double, double>(-1.0, 1.0));
};
template<class System>
void initialize(System system, state_type &inOut, time_type &t, time_type dt)
{
m_coeff.reset();
for(size_t i=0; i<steps+1; ++i)
for(size_t i=0; i<order_value; ++i)
{
do_step(system, inOut, t, dt);
t += dt;
this->do_step(system, inOut, t, dt/static_cast< Time >(order_value));
t += dt/static_cast< Time >(order_value);
};
};
template<class System>
void do_step_impl(System system, const state_type & in, time_type t, state_type & out, time_type &dt, error_storage_type &xerr)
{
size_t eO = m_coeff.m_eo;
m_xerr_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_coeff.predict(t, dt);
if (eO == 1)
{
system(in, m_dxdt.m_v, t);
m_coeff.do_step(m_dxdt.m_v, 1);
}
}
const coeff_type& coeff() const
{
return m_coeff;
out = in;
for(size_t i=0; i<eO; ++i)
{
this->m_algebra.for_each3(out, out, m_coeff.phi[1][i].m_v,
typename Operations::template scale_sum2<double, double>(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<double, double>(1.0, dt*m_coeff.g[eO]));
// error for current order
this->m_algebra.for_each2(xerr[1].m_v, m_coeff.phi[0][eO].m_v,
typename Operations::template scale_sum1<double>(dt*(m_coeff.g[eO]-m_coeff.g[eO-1])));
// error for order below/above
if (eO > 1)
{
this->m_algebra.for_each2(xerr[0].m_v, m_coeff.phi[0][eO-1].m_v,
typename Operations::template scale_sum1<double>(dt*(m_coeff.g[eO-1]-m_coeff.g[eO-2])));
}
this->m_algebra.for_each2(xerr[2].m_v, m_coeff.phi[0][eO+1].m_v,
typename Operations::template scale_sum1<double>(dt*(m_coeff.g[eO+1]-m_coeff.g[eO])));
};
coeff_type& coeff()
{
return m_coeff;
};
const coeff_type& coeff() const {return m_coeff;};
coeff_type & coeff() {return m_coeff;};
wrapped_deriv_type m_dxdt;
void reset() { m_coeff.reset(); };
private:
template< class StateType >
@@ -180,29 +193,35 @@ private:
return adjust_size_by_resizeability( m_dxdt, x, typename is_resizeable<deriv_type>::type() );
};
template< class StateType >
bool resize_xerr_impl( const StateType &x )
{
return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable<deriv_type>::type() );
};
template< class StateType >
bool resize_xnew_impl( const StateType &x )
{
return adjust_size_by_resizeability( m_xnew, x, typename is_resizeable<state_type>::type() );
};
template< class StateType >
bool resize_xerr_impl( const StateType &x )
{
bool resized( false );
aab_type m_aab;
aam_type m_aam;
resizer_type m_dxdt_resizer;
resizer_type m_xerr_resizer;
resizer_type m_xnew_resizer;
wrapped_state_type m_xnew;
wrapped_state_type m_xerr;
for(size_t i=0; i<3; ++i)
{
resized |= adjust_size_by_resizeability( m_xerr[i], x, typename is_resizeable<state_type>::type() );
}
return resized;
};
coeff_type m_coeff;
resizer_type m_dxdt_resizer;
resizer_type m_xnew_resizer;
resizer_type m_xerr_resizer;
wrapped_deriv_type m_dxdt;
wrapped_state_type m_xnew;
error_storage_type m_xerr;
};
}}}
} // odeint
} // numeric
} // boost
#endif
#endif

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

@@ -1,78 +0,0 @@
#ifndef STEPPER_ADAMS_MOULTON_HPP_INCLUDED
#define STEPPER_ADAMS_MOULTON_HPP_INCLUDED
// helper class for controlled_adams_bashforth_moulton
#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/bind.hpp>
#include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
#include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
#include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template<
size_t Steps,
class State,
class Value = double,
class Deriv = State,
class Time = Value,
class Algebra = typename algebra_dispatcher< State >::algebra_type,
class Operations = typename operations_dispatcher< State >::operations_type ,
class Resizer = initially_resizer
>
class adaptive_adams_moulton: public algebra_stepper_base< Algebra , Operations >
{
public:
static const size_t steps = Steps;
typedef unsigned short order_type;
static const order_type order_value = steps + 1;
typedef State state_type;
typedef Value value_type;
typedef Deriv deriv_type;
typedef Time time_type;
typedef Algebra algebra_type;
typedef Operations operations_type;
typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
typedef detail::adaptive_adams_coefficients<order_value - 1, deriv_type, time_type, algebra_type, operations_type> coeff_type;
typedef adaptive_adams_moulton< Steps , State , Value , Deriv , Time , Resizer > stepper_type;
adaptive_adams_moulton( const algebra_type &algebra = algebra_type())
:algebra_stepper_base_type( algebra )
{};
order_type order() const
{
return order_value;
};
void do_step(coeff_type & coeff, const state_type & in,
time_type t, state_type & out, time_type &dt)
{
coeff.poly.add_root(0);
out = in;
// integrating
for(size_t i=0; i<coeff.m_effective_order; ++i)
{
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));
}
};
};
}
}
}
#endif

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

@@ -15,18 +15,59 @@ namespace boost{
namespace numeric{
namespace odeint {
template<
size_t MaxOrder,
class State,
class Algebra = typename algebra_dispatcher< State >::algebra_type
>
class default_order_adjuster
{
public:
typedef State state_type;
typedef state_wrapper< state_type > wrapped_state_type;
typedef Algebra algebra_type;
default_order_adjuster(const algebra_type &algebra = algebra_type())
: m_algebra(algebra)
{};
void adjust_order(size_t &order, const boost::array<wrapped_state_type, 3> & xerr)
{
if(order > 1 && m_algebra.norm_inf(xerr[0].m_v) < 0.5 * m_algebra.norm_inf(xerr[1].m_v))
{
--order;
}
else if (order < MaxOrder && m_algebra.norm_inf(xerr[2].m_v) < 0.5 * m_algebra.norm_inf(xerr[1].m_v))
{
++order;
}
};
private:
algebra_type m_algebra;
};
template<
class ErrorStepper,
class StepAdjuster = detail::pid_step_adjuster<ErrorStepper::order_value,
typename ErrorStepper::state_type,
typename ErrorStepper::time_type,
detail::H211PI>,
detail::H312PID
>,
class OrderAdjuster = default_order_adjuster<ErrorStepper::order_value,
typename ErrorStepper::state_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;
@@ -37,19 +78,34 @@ public:
typedef Resizer resizer_type;
typedef StepAdjuster step_adjuster_type;
typedef OrderAdjuster order_adjuster_type;
typedef controlled_stepper_tag stepper_category;
typedef typename stepper_type::wrapped_state_type wrapped_state_type;
typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
typedef boost::array<wrapped_state_type, 3> error_storage_type;
typedef controlled_adams_bashforth_moulton<ErrorStepper, StepAdjuster, Resizer> controlled_stepper_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_coeff(m_stepper.coeff()),
:m_stepper(),
m_dxdt_resizer(), m_xerr_resizer(), m_xnew_resizer(),
m_step_adjuster(step_adjuster)
m_step_adjuster(step_adjuster), m_order_adjuster()
{};
template<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 System>
controlled_step_result try_step(System system, state_type & inOut, time_type &t, time_type &dt)
{
@@ -67,18 +123,23 @@ public:
controlled_step_result try_step(System system, const state_type & in, time_type &t, state_type & out, time_type &dt)
{
m_xerr_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_stepper.do_step_impl(system, m_coeff, in, t, out, dt, m_xerr.m_v);
m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_stepper.do_step_impl(system, in, t, out, dt, m_xerr);
coeff_type &coeff = m_stepper.coeff();
size_t prevOrder = coeff.m_eo;
m_order_adjuster.adjust_order(coeff.m_eo, m_xerr);
time_type dtPrev = dt;
dt = m_step_adjuster.adjust_stepsize(dt, m_xerr.m_v);
dt = m_step_adjuster.adjust_stepsize(dt, m_xerr[1 + coeff.m_eo - prevOrder].m_v);
if(dt / dtPrev >= 0.9)
if(dt / dtPrev >= step_adjuster_type::threshold())
{
m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
system(out, m_dxdt.m_v, t+dtPrev);
m_coeff.step(m_dxdt.m_v, t+dtPrev);
m_coeff.confirm();
system(out, m_dxdt.m_v, t+dt);
coeff.do_step(m_dxdt.m_v);
coeff.confirm();
t += dtPrev;
return success;
@@ -89,6 +150,8 @@ public:
}
};
void reset() { m_stepper.reset(); };
private:
template< class StateType >
bool resize_dxdt_impl( const StateType &x )
@@ -98,7 +161,13 @@ private:
template< class StateType >
bool resize_xerr_impl( const StateType &x )
{
return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable<state_type>::type() );
bool resized( false );
for(size_t i=0; i<3; ++i)
{
resized |= adjust_size_by_resizeability( m_xerr[i], x, typename is_resizeable<state_type>::type() );
}
return resized;
};
template< class StateType >
bool resize_xnew_impl( const StateType &x )
@@ -107,10 +176,9 @@ private:
};
stepper_type m_stepper;
typename stepper_type::coeff_type &m_coeff;
wrapped_deriv_type m_dxdt;
wrapped_state_type m_xerr;
error_storage_type m_xerr;
wrapped_state_type m_xnew;
resizer_type m_dxdt_resizer;
@@ -118,6 +186,7 @@ private:
resizer_type m_xnew_resizer;
step_adjuster_type m_step_adjuster;
order_adjuster_type m_order_adjuster;
};
}

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

@@ -1,8 +1,7 @@
#ifndef ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED
#define ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_ADAPTIVE_ADAMS_COEFFICIENTS_HPP_INCLUDED
#include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
#include <boost/numeric/odeint/stepper/detail/polynomial.hpp>
#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
@@ -12,8 +11,6 @@
#include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
#include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
#include <boost/numeric/odeint/util/unwrap_reference.hpp>
#include <iostream>
#include <boost/array.hpp>
@@ -25,112 +22,133 @@ namespace detail {
template<
size_t Steps,
class Deriv,
class Time,
class Value = double,
class Time = double,
class Algebra = typename algebra_dispatcher< Deriv >::algebra_type,
class Operations = typename operations_dispatcher< Deriv >::operations_type ,
class Resizer = initially_resizer
>
struct adaptive_adams_coefficients
class adaptive_adams_coefficients
{
public:
static const size_t steps = Steps;
typedef unsigned short order_type;
static const order_type order_value = steps;
typedef Value value_type;
typedef Deriv deriv_type;
typedef Time time_type;
typedef state_wrapper<deriv_type> wrapped_deriv_type;
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 rotating_buffer<wrapped_deriv_type, steps+1> step_storage_type; // +1 for moulton
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, Time, Algebra, Operations, Resizer> aac_type;
typedef detail::Polynomial<steps+2, time_type> poly_type;
typedef adaptive_adams_coefficients<Steps, Deriv, Value, Time, Algebra, Operations, Resizer> aac_type;
adaptive_adams_coefficients(const algebra_type &algebra = algebra_type() )
:poly(), m_effective_order(1), m_resizer(), m_algebra(algebra)
{};
void step(const deriv_type &deriv, const time_type &t)
adaptive_adams_coefficients( const algebra_type &algebra = algebra_type())
:m_eo(1), beta(), phi(), m_ns(0), m_time_storage(),
m_algebra(algebra),
m_phi_resizer()
{
m_resizer.adjust_size( deriv , detail::bind( &aac_type::template resize_tss_impl< deriv_type > , detail::ref( *this ) , detail::_1 ) );
for (size_t i=0; i<order_value+2; ++i)
c[0][i] = 1.0/(i+1);
m_tts[0] = t;
m_tss[0][0].m_v = deriv;
g[0] = c[0][0];
for(size_t i=1; i<m_effective_order; ++i)
beta[0][0] = 1;
beta[1][0] = 1;
};
void predict(time_type t, time_type dt)
{
m_time_storage[0] = t;
if (fabs(m_time_storage[0] - m_time_storage[1] - dt) > 1e-16 || m_eo >= m_ns)
{
time_type dt = t - m_ts[i-1];
m_algebra.for_each3(m_tss[i][0].m_v, m_tss[i-1][0].m_v, m_ss[i-1][0].m_v, typename Operations::template scale_sum2<double, double>(1/dt, -1/dt));
m_ns = 0;
}
else if (m_ns < order_value + 2)
{
m_ns++;
}
for(size_t i=1+m_ns; i<m_eo+1; ++i)
{
beta[0][i] = beta[0][i-1]*(m_time_storage[0] + dt -
m_time_storage[i-1])/(m_time_storage[0] - m_time_storage[i]);
}
for(size_t i=1; i<m_eo+2; ++i)
{
for(size_t j=0; j<m_eo+1; ++j)
{
c[i][j] = c[i-1][j] - c[i-1][j+1]*dt/(m_time_storage[0] + dt - m_time_storage[i-1]);
}
g[i] = c[i][0];
}
};
void do_step(const deriv_type &dxdt, const int o = 0)
{
m_phi_resizer.adjust_size( dxdt , detail::bind( &aac_type::template resize_phi_impl< deriv_type > , detail::ref( *this ) , detail::_1 ) );
phi[o][0].m_v = dxdt;
for(size_t i=1; i<m_eo + 2; ++i)
{
this->m_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<double, double>(1.0, -beta[o][i-1]));
}
};
void confirm()
{
for(size_t i=0; i<steps; ++i)
{
m_ss[i] = m_tss[i];
m_tss[i].rotate();
}
m_ts = m_tts;
m_tts.rotate();
if(m_effective_order < steps+1)
{
++m_effective_order;
}
beta.rotate();
phi.rotate();
m_time_storage.rotate();
};
void reset()
{
poly.reset();
m_effective_order = 1;
m_eo = 1;
};
void pretty_print()
{
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;
}
};
size_t m_eo;
poly_type poly;
time_storage_type m_c;
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;
rotating_buffer<boost::array<value_type, (order_value+1)>, 2> beta; // beta[0] = beta(n)
rotating_buffer<boost::array<wrapped_deriv_type, (order_value + 2)>, 3> phi; // phi[0] = phi(n+1)
boost::array<value_type, order_value + 2> g;
private:
template< class StateType >
bool resize_tss_impl( const StateType &x )
bool resize_phi_impl( const StateType &x )
{
bool resized( false );
for( size_t i=0 ; i<steps+1 ; ++i )
for(size_t i=0; i<(order_value + 2); ++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;
resized |= adjust_size_by_resizeability( phi[0][i], x, typename is_resizeable<deriv_type>::type() );
resized |= adjust_size_by_resizeability( phi[1][i], x, typename is_resizeable<deriv_type>::type() );
resized |= adjust_size_by_resizeability( phi[2][i], x, typename is_resizeable<deriv_type>::type() );
}
return resized;
};
resizer_type m_resizer;
size_t m_ns;
time_storage_type m_time_storage;
boost::array<boost::array<value_type, order_value + 2>, order_value + 2> c;
algebra_type m_algebra;
resizer_type m_phi_resizer;
};
}
@@ -138,4 +156,4 @@ private:
}
}
#endif
#endif

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

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