Как бы я создал пользовательский степпер для интегратора 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.
Спасибо,
Золтан
Попробуй это:
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;
// ...
};
Других решений пока нет …