Hi everyone,
I am trying to solve a system of differential equations using GSL in object oriented programming. I am following the procedure described in:
http://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html
Therefore I need to define functions ('func' and 'jac', see below), the pointers to these functions should then be passed to define the system 'sys'. However, if these functions are defined in classes, they are not of the right type. To solve this, I wrote a wrapper class. This compiles, but still there is something wrong. The functions I am trying to integrate just stay at there initial values. Whereas, when I try exactly the same but not in a class environment, everything works out fine. Any suggestion of how to solve this is welcome!
This is the part of the code that (I think) causes the problems:
template <class T>
class functionWrapper
{
public:
static void initialize(T* funcPtr){_funcPtr = funcPtr;}
static int evaluate_func(double x, const double* y, double* z , void* params)
{
return (*_funcPtr).func(x, y, z, params);
}
static int evaluate_jac(double x, const double* y, double* z, double* w , void* params)
{
return (*_funcPtr).jac(x, y, z, w, params);
}
private:
static T* _funcPtr;
};
template <typename T>
T* functionWrapper<T>::_funcPtr;
class CRGEs
{
...
int func(double, const double*, double*, void*);
int jac(double, const double*, double*, double*, void*);
int integrate();
...
};
int CRGEs::func(double t, const double y[], double f[], void * params)
{
...
}
int CRGEs::jac(double t, const double yuk[], double *dfdyu, double fdtu[],void* params)
{
...
}
int CRGEs::integrate(double y[])
{
...
functionWrapper<CRGEs>::initialize(this);
gsl_odeiv2_system sys = {functionWrapper<CRGEs>::evaluate_func, functionWrapper<CRGEs>::evaluate_jac, dimension, 0};
...
}