Cannot use bulirsch_stoer try_step in place of runge_kutta4 do_step

[iljah]

Compiling the following program gives a lot of errors when using try_step instead of do_step:

#include <boost/numeric/odeint.hpp>
#include <Eigen/Core>
int main() {
    using state_t = std::array<Eigen::Vector3d, 1>;
    //boost::numeric::odeint::runge_kutta4<state_t> stepper;
    boost::numeric::odeint::bulirsch_stoer<state_t> stepper;

    auto tracer = [](const state_t& x, state_t& dxdt, const double /*t*/){
        dxdt[0][0] = dxdt[0][1] = dxdt[0][2] = 1;
    };

    state_t position{{{2, 2, 0}}};
    double time = 0, time_step = 0.01;

    std::cout << position[0].norm() << std::endl;
    //stepper.do_step(tracer, position, time, time_step);
    stepper.try_step(tracer, position, time, time_step);
    std::cout << position[0].norm() << std::endl;

    return 0;
}

The first error looks like:

boost/numeric/odeint/algebra/array_algebra.hpp:283:96: virhe: no matching function for call to ”abs(const value_type&)”
             init = max BOOST_PREVENT_MACRO_SUBSTITUTION ( init , static_cast< result_type >(abs(s[i])) );
                                                                                                ^

If I uncomment runge_kutta4 and do_step everything works fine. The documentation at http://www.boost.org/doc/libs/1_60_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/steppers.html gives no indication that try_step wouldn't work here so is this a bug in odeint or is try_step supposed to be used differently?

[headmyshoulder]

I am not definitely sure, but I think you also need to include

include <boost/numeric/odeint/external/eigen/eigen.hpp>

On 10.03.2016 09:59, Ilja wrote:

Compiling the following program gives a lot of errors when using try_step instead of do_step:

|#include <boost/numeric/odeint.hpp> #include <Eigen/Core> int main() { using state_t = std::array<Eigen::Vector3d, 1>; //boost::numeric::odeint::runge_kutta4 stepper; boost::numeric::odeint::bulirsch_stoer stepper; auto tracer = [](const state_t& x, state_t& dxdt, const double /t/){ dxdt[0][0] = dxdt[0][1] = dxdt[0][2] = 1; }; state_t position{{{2, 2, 0}}}; double time = 0, time_step = 0.01; std::cout << position[0].norm() << std::endl; //stepper.do_step(tracer, position, time, time_step); stepper.try_step(tracer, position, time, time_step); std::cout << position[0].norm() << std::endl; return 0; } |

The first error looks like:

|boost/numeric/odeint/algebra/array_algebra.hpp:283:96: virhe: no matching function for call to ”abs(const value_type&)” init = max BOOST_PREVENT_MACRO_SUBSTITUTION ( init , static_cast< result_type

(abs(s[i])) ); ^ |

If I uncomment runge_kutta4 and do_step everything works fine. The documentation at http://www.boost.org/doc/libs/1_60_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/steppers.html gives no indication that try_step wouldn't work here so is this a bug in odeint or is try_step supposed to be used differently?

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/191.

[iljah]

Now I get a much shorter error:

boost/numeric/odeint/algebra/array_algebra.hpp:283:57: error: invalid static_cast from type 'const Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs_op<double>, const Eigen::Matrix<double, 3, 1> >' to type 'result_type {aka double}'
             init = max BOOST_PREVENT_MACRO_SUBSTITUTION ( init , static_cast< result_type >(abs(s[i])) );
                                                         ^

while runge_kutta4 version still works.

[headmyshoulder]

Ok, it seems that an overload for abs is missing. I will have a look today at this problem and come back with a solution.

On 10.03.2016 10:23, Ilja wrote:

Now I get a much shorter error:

|boost/numeric/odeint/algebra/array_algebra.hpp:283:57: error: invalid static_cast from type 'const Eigen::CwiseUnaryOpEigen::internal::scalar_abs_op<double, const Eigen::Matrix<double, 3, 1> >' to type 'result_type {aka double}' init = max BOOST_PREVENT_MACRO_SUBSTITUTION ( init , static_cast< result_type

(abs(s[i])) ); ^ |

while runge_kutta4 version still works.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/191#issuecomment-194749131.

[headmyshoulder]

Ok, this issue is a bit more complicated. The problem is, that the existing overload of abs() is expected to calculate the absolute value component-wise which is not was is needed here if I see this correctly.

It might be possible to get this working by providing the correct overloads for several classes. I will try to find a solution. Btw. do you need to use an array of such a vector? Maybe it is easier to solve your problem using vector< double> or something similar.

On 10.03.2016 10:23, Ilja wrote:

Now I get a much shorter error:

|boost/numeric/odeint/algebra/array_algebra.hpp:283:57: error: invalid static_cast from type 'const Eigen::CwiseUnaryOpEigen::internal::scalar_abs_op<double, const Eigen::Matrix<double, 3, 1> >' to type 'result_type {aka double}' init = max BOOST_PREVENT_MACRO_SUBSTITUTION ( init , static_cast< result_type

(abs(s[i])) ); ^ |

while runge_kutta4 version still works.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/191#issuecomment-194749131.

[mariomulansky]

Indeed, using a std::array<Eigen::Vector3d, 1> as state type is currently not working due to some contradictory uses of abs. In your specific case, you should be able to just use Eigen::Vector3d? However, I see some use in being able to work with arrays of Eigen vectors so we should try to fix this.

[iljah]

My use case is tracing a vector field so the state has to be an array of some sort I think.

A long time ago I noticed that I couldn't use types which don't have [] operator as a state, i.e. a double, but that seems to have been fixed now. In this case swithing state_t from array to just vector3d works, thanks!

[mariomulansky]

It was always possible to use double or Eigen vectors as state types, but in earlier version you needed to select the correct algebra backend yourself. At some point we introduced some automatic mechanism that's why it now just works out of the box :)

However, your original problem is still valid and should be addressed, we will look into that.