Skip to content

Commit

Permalink
Improvement to Order Selection, Error approximation in the ABM stepper (
Browse files Browse the repository at this point in the history 
#218)

* improves the order selection and modifies the error estimation accordingly

- assumes constant stepsize for the next step to approximate error
- moves the complete order adjustment to the class order_adjustment
- slight changes to adaptive_adams_coefficients

* changed the commit according to the requests and comments
  • Loading branch information
@mariomulansky
Valentin Hartmann authored and mariomulansky committed Sep 25, 2017
1 parent af59be4 commit 540f46f
Show file tree
Hide file tree
Showing 3 changed files with 126 additions and 91 deletions.
43 changes: 10 additions & 33 deletions include/boost/numeric/odeint/stepper/adaptive_adams_bashforth_moulton.hpp
View file
Expand Up @@ -49,7 +49,6 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O

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 typename algebra_stepper_base_type::algebra_type algebra_type;
Expand All @@ -74,31 +73,29 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O
{
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);
do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr.m_v);
boost::numeric::odeint::copy( m_xnew.m_v , inOut);
};

template< class System >
void do_step(System system, const state_type &in, time_type t, state_type &out, time_type dt )
{
do_step(system, in, t, out, dt, m_xerr);
do_step(system, in, t, out, dt, m_xerr.m_v);
};

template< class System >
void do_step(System system, state_type &inOut, time_type t, time_type dt, state_type &xerr)
{
m_xnew_resizer.adjust_size( inOut , detail::bind( &stepper_type::template resize_xnew_impl< state_type > , detail::ref( *this ) , detail::_1 ) );

do_step(system, inOut, t, m_xnew.m_v, dt, m_xerr[1].m_v);
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_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);
do_step_impl(system, in, t, out, dt, xerr);

system(out, m_dxdt.m_v, t+dt);
m_coeff.do_step(m_dxdt.m_v);
Expand Down Expand Up @@ -142,15 +139,15 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O
};

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)
void do_step_impl(System system, const state_type & in, time_type t, state_type & out, time_type &dt, state_type &xerr)
{
size_t eO = m_coeff.m_eo;

m_xerr_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_dxdt_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );

m_coeff.predict(t, dt);
if (eO == 1)
if (m_coeff.m_steps_init == 1)
{
system(in, m_dxdt.m_v, t);
m_coeff.do_step(m_dxdt.m_v, 1);
Expand All @@ -170,22 +167,8 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O
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
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])));
}

// error for order above
if(m_coeff.m_steps_init > eO)
{
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])));
}
this->m_algebra.for_each2(xerr, m_coeff.phi[0][eO].m_v,
typename Operations::template scale_sum1<double>(dt*(m_coeff.g[eO])));
};

const coeff_type& coeff() const { return m_coeff; };
Expand All @@ -208,13 +191,7 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O
template< class StateType >
bool resize_xerr_impl( const StateType &x )
{
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;
return adjust_size_by_resizeability( m_xerr, x, typename is_resizeable<state_type>::type() );
};

coeff_type m_coeff;
Expand All @@ -225,7 +202,7 @@ class adaptive_adams_bashforth_moulton: public algebra_stepper_base< Algebra , O

wrapped_deriv_type m_dxdt;
wrapped_state_type m_xnew;
error_storage_type m_xerr;
wrapped_state_type m_xerr;
};

} // odeint
Expand Down
116 changes: 74 additions & 42 deletions include/boost/numeric/odeint/stepper/controlled_adams_bashforth_moulton.hpp
View file
Expand Up @@ -39,41 +39,58 @@ class default_order_adjuster
: m_algebra( algebra )
{};

void adjust_order(size_t &order, const boost::array<wrapped_state_type, 3> & xerr)
size_t adjust_order(size_t order, size_t init, boost::array<wrapped_state_type, 4> &xerr)
{
using std::abs;

value_type errm = abs(m_algebra.norm_inf(xerr[0].m_v));
value_type errc = abs(m_algebra.norm_inf(xerr[1].m_v));
value_type errp = abs(m_algebra.norm_inf(xerr[2].m_v));
value_type errc = abs(m_algebra.norm_inf(xerr[2].m_v));

if(order == 1)
{
if(errp < errc)
{
order ++;
}
value_type errm1 = 3*errc;
value_type errm2 = 3*errc;

if(order > 2)
{
errm2 = abs(m_algebra.norm_inf(xerr[0].m_v));
}
else if(order == MaxOrder)
if(order >= 2)
{
if(errm < errc)
{
order --;
}
errm1 = abs(m_algebra.norm_inf(xerr[1].m_v));
}
else

size_t o_new = order;

if(order == 2 && errm1 <= 0.5*errc)
{
o_new = order - 1;
}
else if(order > 2 && errm2 < errc && errm1 < errc)
{
o_new = order - 1;
}

if(init < order)
{
return order+1;
}
else if(o_new == order - 1)
{
return order-1;
}
else if(order <= MaxOrder)
{
if(errm < errc &&
errm < errp)
value_type errp = abs(m_algebra.norm_inf(xerr[3].m_v));

if(order > 1 && errm1 < errc && errp)
{
order--;
return order-1;
}
else if(errp < errm &&
errp < errc)
else if(order < MaxOrder && errp < (0.5-0.25*order/MaxOrder) * errc)
{
order++;
return order+1;
}
}

return order;
};
private:
algebra_type m_algebra;
Expand Down Expand Up @@ -118,7 +135,7 @@ class controlled_adams_bashforth_moulton

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 boost::array< wrapped_state_type , 4 > error_storage_type;

typedef typename stepper_type::coeff_type coeff_type;
typedef controlled_adams_bashforth_moulton< ErrorStepper , StepAdjuster , OrderAdjuster , Resizer > controlled_stepper_type;
Expand Down Expand Up @@ -162,39 +179,54 @@ class controlled_adams_bashforth_moulton
m_xerr_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_xerr_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_stepper_type::template resize_dxdt_impl< state_type > , detail::ref( *this ) , detail::_1 ) );

m_stepper.do_step_impl(system, in, t, out, dt, m_xerr);
m_stepper.do_step_impl(system, in, t, out, dt, m_xerr[2].m_v);

coeff_type &coeff = m_stepper.coeff();

time_type dtPrev = dt;
size_t prevOrder = coeff.m_eo;

if(coeff.m_steps_init > coeff.m_eo)
{
m_order_adjuster.adjust_order(coeff.m_eo, m_xerr);
dt = m_step_adjuster.adjust_stepsize(coeff.m_eo, dt, m_xerr[1 + coeff.m_eo - prevOrder].m_v, in, m_stepper.dxdt());
}
else
{
if(coeff.m_eo < order_value)
{
coeff.m_eo ++;
}
dt = m_step_adjuster.adjust_stepsize(coeff.m_eo, dt, m_xerr[1].m_v, in, m_stepper.dxdt());
}
dt = m_step_adjuster.adjust_stepsize(coeff.m_eo, dt, m_xerr[2].m_v, out, m_stepper.dxdt() );

if(dt / dtPrev >= step_adjuster_type::threshold())
{
system(out, m_dxdt.m_v, t+dt);
system(out, m_dxdt.m_v, t+dtPrev);

coeff.do_step(m_dxdt.m_v);
coeff.confirm();

t += dtPrev;

size_t eo = coeff.m_eo;

// estimate errors for next step
double factor = 1;
algebra_type m_algebra;

m_algebra.for_each2(m_xerr[2].m_v, coeff.phi[1][eo].m_v,
typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo])));

if(eo > 1)
{
m_algebra.for_each2(m_xerr[1].m_v, coeff.phi[1][eo-1].m_v,
typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo-1])));
}
if(eo > 2)
{
m_algebra.for_each2(m_xerr[0].m_v, coeff.phi[1][eo-2].m_v,
typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo-2])));
}
if(eo < order_value && coeff.m_eo < coeff.m_steps_init-1)
{
m_algebra.for_each2(m_xerr[3].m_v, coeff.phi[1][eo+1].m_v,
typename operations_type::template scale_sum1<double>(factor*dt*(coeff.gs[eo+1])));
}

// adjust order
coeff.m_eo = m_order_adjuster.adjust_order(coeff.m_eo, coeff.m_steps_init-1, m_xerr);

return success;
}
else
{
coeff.m_eo = prevOrder;
return fail;
}
};
Expand All @@ -212,7 +244,7 @@ class controlled_adams_bashforth_moulton
{
bool resized( false );

for(size_t i=0; i<3; ++i)
for(size_t i=0; i<m_xerr.size(); ++i)
{
resized |= adjust_size_by_resizeability( m_xerr[i], x, typename is_resizeable<state_type>::type() );
}
Expand Down

0 comments on commit 540f46f

Please sign in to comment.