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

C++ class member function to GSL ODE solver

本文关键字：ODE GSL 成员 函数 C++      更新时间：2023-10-16

显然，我正在努力解决一个众所周知的问题。我有一个由类成员函数定义的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%.3ft%.3ft%.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，会导致所有的分段错误(但我认为我在实现中做错了什么，因为那篇文章中的指示非常模糊)。如果能为这个问题提供一个一劳永逸的、尽可能简单的工作解决方案(可能不使用boost库)，那就太好了。

类似这样的东西：

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

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

相关文章：