gsl_odeiv2在c + +类:模板包装int(...)ode
我目前在我的类中使用gsl_odeiv2方法来解微分方程。但是由于众所周知的成员函数问题,我无法在类中定义我的ode-system。我目前使用的一种变通方法: 我在一个全局命名空间定义我的颂歌:gsl_odeiv2在c + +类:模板包装int(...)ode
ODE.hpp:
#include "EoS.hpp"
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
namespace ODEs
{
struct tov_eq_params {EoS *eos;};
int tov_eq(double, const double *, double *, void *);
}
ODE.cpp:
namespace ODEs {
int tov_eq(double r, const double *PM, double *dPdM, void *p) {
struct tov_eq_params * params = (struct tov_eq_params *)p;
EoS *eos = (params->eos);
...
return GSL_SUCCESS
}
}
,并使用指针coustom类型(类EOS)作为参数的对象。内部类解决我的颂歌我使用:
...
struct tov_eq_params comp_tov_params = {(this->star_eos)};
gsl_odeiv2_evolve *comp_tov_evolve = gsl_odeiv2_evolve_alloc(3);
gsl_odeiv2_system comp_tov_system = {tov_eq, NULL, 3,&comp_tov_params};
...
initalise我的系统。这工作正常,但有点混乱,因为我需要在全局命名空间中声明我的微分方程。
我知道有可能使用gsl_functions的模板包装stackoverflow.com/questions/.../how-to-avoid-static-member-function-when-using-gsl-with-c/...在C++类中使用它们。我实际上使用了那里描述的包装来为我的类中的gsl_integration方法定义函数,它完美地工作,并且更简洁,代码更少。例如:我可以用我的star_eos从上面direcly对象在函数内部:
auto dBf = [=](double r)->double{
return 4 * M_PI * gsl_pow_2(r) * (this->star_eos)->nbar(this->P(r)) * sqrt(this->expLambda(r))* 1e54;
};
gsl_function_pp<decltype(dBf)> dBfp(dBf);
gsl_function *dB = static_cast<gsl_function*>(&dBfp);
我试着写这样的INT(双R,常量双* PM,模板包装双* DPDM,无效* P )函数gsl_odeiv2_system需要,但我失败了,因为我是C++的新手,并没有完全理解它的模板/ static_cast机制。
是否有人将gsl_odeiv方法及其ode系统与模板包装器一起使用?或者可以有人想出一个模板类似于上面描述的gsl_functions但是int(...)ode。
思考我如何在全局命名空间中设置微分方程时使用它,我找到了解决方案。我现在有一个工作包装。在一个全局命名空间,我有以下:
//gsl_wrapper.hpp
#include <iostream>
#include <vector>
#include <functional>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
namespace gsl_wrapper {
class ode_System{
public:
ode_System(int);
int dim;
std::function<double (double, const double *, double *, int)> *df;
};
struct ode_struct {ode_System *ode;};
int ode(double, const double *, double *, void *);
}
//gsl_wrapper.cpp
#include "../../include/gsl_wrapper.hpp"
namespace gsl_wrapper {
ode_System::ode_System(int dim) {
this->dim=dim;
}
int ode(double r, const double *f, double *df, void *p) {
struct ode_struct * params = (struct ode_struct *)p;
ode_System *odeFunc = (params->ode);
int dim = odeFunc->dim;
std::function<double (double, const double *, double *, int)> dfeq=*(odeFunc->df);
for(int i=0;i<dim;i++){
df[i] = dfeq(r,f,df,i);
}
return GSL_SUCCESS;
}
};
所以我bassically有存储在我的新类ode_System我所有的具体信息,其中有一个INT昏暗到指定的系统尺寸和指针,以便一std :: function对象。这个对象代表了数学微分方程组。
里面我的类,其中,我要解决使用gsl_odeiv2一微分方程,我定义系统使用lambda函数:
std::function<double (double, const double *, double *, int)> dPMeq = [=](double r , const double * PM, double *dPdM, int i)->double{
double df;
switch (i)
{
case 0:
df = ... // df_1/dr
break;
case 1:
df = ... // df_2/dr
break;
case 2:
df = ... // df_3/dr
break;
default:
GSL_ERROR_VAL ("comp_tov_eq:", GSL_FAILURE,GSL_FAILURE);
df = 0;
}
return df;
};
上述系统代表3次微分方程的耦合系统。然后我声明一个ode_System对象具有正确的尺寸,并将其功能指针df设置为我定义的系统。然后,我只需要一个结构具有参考该系统完成:我可以用我的类中定义我微分方程gsl_odeiv2_system:
ode_System tov(3);
tov.df= &dPMeq;
struct ode_struct comp_tov_params = {&tov};
gsl_odeiv2_evolve *comp_tov_evolve = gsl_odeiv2_evolve_alloc(3);
gsl_odeiv2_system comp_tov_system = {ode, NULL, 3, &comp_tov_params};
...
至于我可以告诉这个作品一样好(或坏)作为我在我的问题中提出的实现。它可以使用一些清理,但原则上这对我来说工作得很好。
但是,如果有人知道更好的方式来做到这一点,请随时分享!