GSL ODE 求解器的 C++ 类成员函数

Question

显然，我正在努力解决一个众所周知的问题。我有一个由类成员函数定义的 ODE 系统，我想通过其中一个 GSL 求解器对其进行求解/积分。即，假设我的模型定义为:

class my_model{
...
public:
    int NEQ = 4;
    double y[4], dydt[4];
    double params[25]; 
    int ode_gsl(double t, const double y[], double dydt[], void * params);
    ...
};

int my_model::int ode_gsl(double t, const double y[], double dydt[], void * params){
    dydt[0] = params[1]*params[0]*y[0] - y[1];
    dydt[1] = ...;
    ...
    return GSL_SUCCESS;    
}

然后在我的集成例程中:

int main(){
    my_model chemosc;

   // Parameters for the Integrator
   double HSTART = 1e-3;
   double ATOL = 1e-8;
   double RTOL = 1e-6;

   // Instantiate GSL solver
   gsl_odeiv2_system sys = {&chemosc.ode_gsl, nullptr, chemosc.NEQ, chemosc.params};
   gsl_odeiv2_driver * d;
   d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, HSTART, ATOL, RTOL);

   double twin[2] = {0.,60.};
   double dt = 1e-3;
   double t = twin[0], t1 = twin[1];
   long int NSTEP = (long int)((t1-t)/dt)+1; // +1 if you start counting from zero...
   int NEQ = 4;
   long int NUMEL = (NEQ+1)*NSTEP; // number of elements for solution

   int i = 0,j;
   do{
      double ti = t + dt; // t is automatically updated by the driver
      printf("\n%.3f\t%.3f\t%.3f t%.3f",astro.y[0],astro.y[1],astro.y[2],astro.y[3]);
      int status = gsl_odeiv2_driver_apply (d, &t, ti, astro.y);
      ...
      }
...
}

编译上述代码会产生一个众所周知的错误，即 GSL 需要指向函数的指针，而我传递的是指向成员函数的指针，即:

error: cannot convert ‘int (chemosc::*)(double, const double*, double*, void*)’ to ‘int (*)(double, const double*, double*, void*)’

我发现了几个与此主题相关的问题:Q1 , Q2 , Q3 , Q4 , Q5 , Q6 ，但严重的是，答案很难遵循。将我的成员函数声明为 static 有一个缺点，即编译器要求我也将所有成员参数声明为 static。按照建议使用 static_cast here ，导致所有的段错误(但我假设我在实现中做错了什么，因为那篇文章中的方向非常神秘)。对于这个问题(可能不使用 boost 库)有一个一劳永逸的清晰且尽可能简单的工作解决方案会很好。

Answer 1

像这样:

class my_model
{
    private: int
    ode_gsl_impl(double const t, double const * const y, double * const dydt);

    public: static int
    ode_gsl(double const t, double const * const y, double * const dydt, void * const opaque)
    {
        assert(opaque);
        return(static_cast<my_model *>(opaque)->ode_gsl_impl(t, y, dydt));
    }
};

gsl_odeiv2_system sys
{
    &my_model::ode_gsl
,   nullptr
,   chemosc.NEQ
,   reinterpret_cast<void *>(::std::addressof(chemosc))
};

我想提一下，您的原始代码在回调参数名称和类字段名称之间存在名称冲突。