Define custom stepper for dynamic array in boost

Question

How would I whip up a custom stepper for the ODE integrator in boost? I know how to do that for an array whose size is known at compile time. A simple implementation is like this

#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;

typedef boost::array<double, 2> state_type;

template <size_t N>
class euler_stepper {
public:
    typedef double value_type;
    typedef double time_type;
    typedef unsigned short order_type;

    typedef boost::numeric::odeint::stepper_tag stepper_category;

    static order_type order(void) { return(1); }

    template<class System>
    void do_step(System system, state_type &x, time_type t, time_type dt) const {
        state_type der;
        system(x , der);
        for(size_t i=0 ; i<x.size(); i++) {
            x[i] += dt*der[i];
        }
    }
};

struct rhs {
    void operator()(const state_type &x, state_type &dxdt) const {
        for(int i=0; i < x.size(); i++) dxdt[i] = -x[i];
    }
};

struct observer {
    void operator()(const state_type &x , double t) const {
        for(int i=0; i < x.size(); i++) cout << x[i] << '\t';
        cout << endl;
    }
};

int main(int argc , char **argv) {  
    double dt = 0.01;
    state_type x = {{2.0, 1.0}};
    integrate_const(euler_stepper<2>(), rhs(), x, 0.0, 0.1, dt, observer());
    return 0;
}

But how should I change the implementation, if I wanted to work with a state_type like this

typedef vector<complex<double>> state_type;

Basically, I would like to pass the size of the state vector as an argument to main.

Thanks,

Zoltán

There are several similar examples in the odeint documentation. Are these examples sufficient for you? Usually, there is no need to implement the steppers yourself. I don't see any special requirements why you need to do it here. — headmyshoulder, Nov 06 '15 at 22:24
I would like to solve a stochastic equation, and there are no ready-to-use steppers for that. The example in the documentation is based on arrays, not vectors. — v923z, Nov 07 '15 at 06:35
Ok, I see. But look for the ode stepper. They are usually parametrized with the state type, which can be either a vector or an array. And it is easy to do the same for your example. — headmyshoulder, Nov 07 '15 at 09:30

score 2 · Accepted Answer · answered Nov 07 '15 at 09:34

Try this:

template < typename State >
class euler_stepper {
public:
    using state_type = State;
    using value_type = typename state_type::value_type;
    using time_type = value_type;
    using order_type = unsigned short;

    using stepper_category = boost::numeric::odeint::stepper_tag;

    static order_type order(void) { return(1); }

    template<class System>
    void do_step( System system , state_type &x , time_type t , time_type dt ) const
    {
        state_type der;
        system(x , der);
        for(size_t i=0 ; i<x.size(); i++) {
            x[i] += dt*der[i];
        }
    }
};

You can also fix the state type for your solver, for example:

class euler_stepper
{
public:
    using state_type = vector< complex< double > >;
    using value_type = double;
    using time_type = double;
    // ...
};

Many thanks, this indeed fixes the problem. The only thing missing was that one has to reserve space for the derivative, otherwise, the code segfaults. But that was a simple ```state_type der(x.size());``` — v923z, Nov 08 '15 at 08:25

Define custom stepper for dynamic array in boost

1 Answers1