2

我将如何在 boost 中为 ODE 积分器创建一个自定义步进器？我知道如何对在编译时知道大小的数组执行此操作。一个简单的实现是这样的

#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;
}

但是，如果我想使用这样的 state_type，我应该如何更改实现

typedef vector<complex<double>> state_type;

基本上，我想将状态向量的大小作为参数传递给 main。

谢谢，

佐尔坦

4

1 回答 1

2

试试这个：

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];
        }
    }
};

您还可以修复求解器的状态类型，例如：

class euler_stepper
{
public:
    using state_type = vector< complex< double > >;
    using value_type = double;
    using time_type = double;
    // ...
};
于 2015-11-07T09:34:06.780 回答