Search code examples
c++boostodeint

boost odeint: controlled stepper with custom class and vector space algebra


I am using boost odeint to evolve a differential equation. The time and the value are given by a wrapper class around a plain number (doing Gaussian error propagation).

using Time = Number<double>;
using Value = Number<double>;

The state and its derivative are given by a custom class with the declaration

class State : boost::addable<State>, boost::multipliable<State, Time>, boost::dividable<State>
{
public:
    void Add(Group const& group, Value alpha);
    Value Get(Group const& group) const;
    State& operator*=(Time const& time);
    State& operator+=(State const& state);
    State& operator/=(State const& state);
    State abs() const;
    Value norm_inf() const;
private:
    std::map<Group, Value> map_;
};

State abs(State const& state);

namespace boost
{
namespace numeric
{
namespace odeint
{
template<>
struct boost::numeric::odeint::vector_space_norm_inf<State> {
    using result_type = Value;
    result_type operator()(State const& state) const;
};
}
}
}

using Derivative = State;

I use the Runge-Kutta-Dopri-5 Stepper with the Vector space algebra

using Stepper = boost::numeric::odeint::runge_kutta_dopri5<State, Value, Derivative, Time, boost::numeric::odeint::vector_space_algebra>;

And I use the whole think like

auto start_state = GetStart();
auto system = GetSystem();
auto min = GetMin();
auto max = GetMax();
auto dt = GetDt;
auto observer = GetObserver();
auto stepper = boost::numeric::odeint::make_controlled<Stepper>(1e-6, 1e-6);
boost::numeric::odeint::integrate_adaptive(Stepper(), system, start_state, min, max, dt, observer);

This works like a charm. But as you might notice I do not use the controlled stepper I define in the second last line, but I use the default stepper. The moment I try to use the stepper variable instead of a default Stepper() in the last line I get follwoing compilation error:

In file included from /home/me/Integrate.cpp:2:
In file included from /usr/include/boost/numeric/odeint.hpp:35:
/usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:85:17: error: no matching constructor for initialization of 'typename operations_type::rel_error<value_type>' (aka 'rel_error<me::Number<double> >')
                typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * abs(get_unit_value( dt )) ) );
                ^                                                           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:768:50: note: in instantiation of function template specialization 'odeint::default_error_checker<me::Number<double>, odeint::vector_space_algebra, odeint::default_operations>::error<me::State, me::State, me::State, me::Number<double> >' requested here
        value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt_in , m_xerr.m_v , dt );
                                                 ^
/usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:716:38: note: in instantiation of function template specialization 'odeint::controlled_runge_kutta<odeint::runge_kutta_dopri5<me::State, me::Number<double>, me::State, me::Number<double>, odeint::vector_space_algebra, odeint::default_operations, odeint::initially_resizer>, odeint::default_error_checker<me::Number<double>, odeint::vector_space_algebra, odeint::default_operations>, odeint::default_step_adjuster<me::Number<double>, me::Number<double> >, odeint::initially_resizer, odeint::explicit_error_stepper_fsal_tag>::try_step<(lambda at /home/me/Integrate.cpp:28:19), me::State, me::State, me::State, me::State>' requested here
        controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , m_dxdtnew.m_v , dt );
                                     ^
/usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:899:16: note: in instantiation of function template specialization 'odeint::controlled_runge_kutta<odeint::runge_kutta_dopri5<me::State, me::Number<double>, me::State, me::Number<double>, odeint::vector_space_algebra, odeint::default_operations, odeint::initially_resizer>, odeint::default_error_checker<me::Number<double>, odeint::vector_space_algebra, odeint::default_operations>, odeint::default_step_adjuster<me::Number<double>, me::Number<double> >, odeint::initially_resizer, odeint::explicit_error_stepper_fsal_tag>::try_step<(lambda at /home/me/Integrate.cpp:28:19), me::State, me::State>' requested here
        return try_step( system , x , m_dxdt.m_v , t , dt );
               ^
/usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:617:16: note: in instantiation of function template specialization 'odeint::controlled_runge_kutta<odeint::runge_kutta_dopri5<me::State, me::Number<double>, me::State, me::Number<double>, odeint::vector_space_algebra, odeint::default_operations, odeint::initially_resizer>, odeint::default_error_checker<me::Number<double>, odeint::vector_space_algebra, odeint::default_operations>, odeint::default_step_adjuster<me::Number<double>, me::Number<double> >, odeint::initially_resizer, odeint::explicit_error_stepper_fsal_tag>::try_step_v1<(lambda at /home/me/Integrate.cpp:28:19), me::State>' requested here
        return try_step_v1( system , x , t , dt );
               ^
/usr/include/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp:103:22: note: in instantiation of function template specialization 'odeint::controlled_runge_kutta<odeint::runge_kutta_dopri5<me::State, me::Number<double>, me::State, me::Number<double>, odeint::vector_space_algebra, odeint::default_operations, odeint::initially_resizer>, odeint::default_error_checker<me::Number<double>, odeint::vector_space_algebra, odeint::default_operations>, odeint::default_step_adjuster<me::Number<double>, me::Number<double> >, odeint::initially_resizer, odeint::explicit_error_stepper_fsal_tag>::try_step<(lambda at /home/me/Integrate.cpp:28:19), me::State>' requested here
            res = st.try_step( system , start_state , start_time , dt );
                     ^
/usr/include/boost/numeric/odeint/integrate/integrate_adaptive.hpp:42:20: note: in instantiation of function template specialization 'odeint::detail::integrate_adaptive<odeint::controlled_runge_kutta<odeint::runge_kutta_dopri5<me::State, me::Number<double>, me::State, me::Number<double>, odeint::vector_space_algebra, odeint::default_operations, odeint::initially_resizer>, odeint::default_error_checker<me::Number<double>, odeint::vector_space_algebra, odeint::default_operations>, odeint::default_step_adjuster<me::Number<double>, me::Number<double> >, odeint::initially_resizer, odeint::explicit_error_stepper_fsal_tag>, (lambda at /home/me/Integrate.cpp:28:19), me::State, me::Number<double>, (lambda at /home/me/Integrate.cpp:48:21)>' requested here
    return detail::integrate_adaptive(
                   ^
/usr/include/boost/numeric/odeint/algebra/default_operations.hpp:435:9: note: candidate constructor not viable: no known conversion from 'me::State' to 'me::Number<double>' for 4th argument
        rel_error( Fac1 eps_abs , Fac1 eps_rel , Fac1 a_x , Fac1 a_dxdt )
        ^
/usr/include/boost/numeric/odeint/algebra/default_operations.hpp:431:12: note: candidate constructor (the implicit move constructor) not viable: requires 1 argument, but 4 were provided
    struct rel_error
           ^
/usr/include/boost/numeric/odeint/algebra/default_operations.hpp:431:12: note: candidate constructor (the implicit copy constructor) not viable: requires 1 argument, but 4 were provided

I can see that odeint::rel_error<me::Number<double>> is supposed to be of the Number type, but the product m_a_dxdt * abs(get_unit_value( dt )), which is used as the 4th parameter of the rel_error constructor is of State type. I assume that I have made a mistake in implementing the vector space algebra for the State class, but I can not find where I have done the mistake.

If I replace the the Time and Value types with plain double's, the error checker works again.


Solution

  • It seems like the instruction of odeint contains a mistake. If the return type of the abs function is Value instead of State as one would naively expect it to be the compilation error vanishes.