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

【问题标题】：C++ class member function to GSL ODE solverGSL ODE 求解器的 C++ 类成员函数
【发布时间】：2018-04-13 12:48:22
【问题描述】：

我正在努力解决 - 显然 - 一个众所周知的问题。我有一个由类成员函数定义的 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_cast 使用here，会导致所有的分段错误（但我认为我在实现中做错了什么，因为该帖子中的方向非常神秘）。对于这个问题，最好有一个一劳永逸的、尽可能简单的解决方案（可能不使用boost库）。

【问题讨论】：

函数签名似乎包含void * 参数，您可以使用该参数将指针传递给类实例（以及您需要的任何参数）。因此，您可以声明一个静态成员函数（或非成员函数），然后将 void * 参数转换为类实例，然后调用将执行所有计算的类成员函数。
@VTT 是的，我有点明白这是继续的方式。但是你能给我来示例代码吗？ stackoverflow.com/questions/10593726/… 中提到了您所说的内容，但同样，它与积分器有关，而不是与 ODE 求解器有关。而且解释仍然有些混乱。

标签： c++ class pointers ode gsl

【解决方案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))
};

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

【讨论】：

是的。它就像一个魅力。谢谢。非常优雅的技巧，非常简洁易懂。