GSL ODE求解器的C ++类成员函数 [英] C++ class member function to GSL ODE solver

问题描述

我正在努力-显然-一个众所周知的问题.我有一个由类成员函数定义的ODE系统，我想通过一个GSL求解器对其进行求解/集成.也就是说，假设我的模型是通过以下方式定义的:

I am struggling with -- apparently -- a well known problem. I have an ODE system which is defined by a class member function, and I want to solve/integrate it by one of the GSL solvers. I.e., say my model is defined by:

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

然后在我的集成例程中:

Then in my integration routine:

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需要一个指向函数的指针，而我正在传递一个指向成员函数的指针，即:

Compiling the above code creates the somehow well known error that GSL requires a pointer-to-function whereas I am passing a pointer-to-member function, i.e.:

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

我发现了与此主题相关的几个问题: Q1 ，第三季度， Q5 ，这里，会导致所有细分错误(但我认为我在实现中做错了，因为那篇文章中的说明很含糊).对于这个问题，有一个一目了然的清晰且尽可能简单的工作解决方案(可能无需使用boost库)，将是很好的选择.

I found several questions related to this topic: Q1, Q2, Q3, Q4, Q5, Q6, but seriously, the answer there are hard to follow. Declaring as static my member function, has the downside that the compiler requires me to also declare as static all the member parameters. Using static_cast as suggested here, result in all bunch of segmentation errors (but I assume that I did something wrong in the implementation, inasmuch as the directions in that post are quite cryptic). It would be nice to have a once-for-all clear and as simple as possible working solution for this problem (possibly without using boost libraries).

推荐答案

类似以下内容:

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

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

I would like to mention that your original code suffers from name collisions between callback parameter names and class field names.

这篇关于GSL ODE求解器的C ++类成员函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！