Buscador desde hace mucho tiempo, primera vez que pregunta aquí. He escrito una serie de scripts para realizar varios métodos de integración numérica 1D y los he compilado en una biblioteca. Me gustaría que esa biblioteca sea lo más flexible posible con respecto a lo que es capaz de integrar.

Aquí incluyo un ejemplo: un ejemplo de regla trapezoidal muy simple donde paso un puntero a la función que se va a integrar.

// Numerically integrate (*f) from a to b
// using the trapezoidal rule.
double trap(double (*f)(double), double a, double b) {
  int N = 10000;
  double step = (b-a)/N;
  double s = 0;
  for (int i=0; i<=N; i++) {
    double xi = a + i*step;
    if (i == 0 || i == N) { s += (*f)(xi); }
    else { s += 2*(*f)(xi); }
  }
  s *= (b-a)/(2*N);
  return s;
}

Esto funciona muy bien para funciones simples que solo toman un argumento. Ejemplo:

double a = trap(sin,0,1);

Sin embargo, a veces es posible que desee integrar algo que tenga más parámetros, como un polinomio cuadrático. En este ejemplo, los coeficientes serían definidos por el usuario antes de la integración. Código de ejemplo:

// arbitrary quadratic polynomial
double quad(double A, double B, double C, double x) {
  return (A*pow(x,2) + B*x + C);
}

Idealmente, podría hacer algo como esto para integrarlo:

double b = trap(quad(1,2,3),0,1);

Pero claramente eso no funciona. He solucionado este problema definiendo una clase que tiene los coeficientes como miembros y la función de interés como una función miembro:

class Model {
  double A,B,C;
public:
  Model() { A = 0; B = 0; C = 0; }
  Model(double x, double y, double z) { A = x; B = y; C = z; }
  double func(double x) { return (A*pow(x,2)+B*x+C); }
};

Sin embargo, entonces mi función de integración necesita cambiar para tomar un objeto como entrada en lugar de un puntero a función:

// Numerically integrate model.func from a to b
// using the trapezoidal rule.
double trap(Model poly, double a, double b) {
  int N = 10000;
  double step = (b-a)/N;
  double s = 0;
  for (int i=0; i<=N; i++) {
    double xi = a + i*step;
    if (i == 0 || i == N) { s += poly.func(xi); }
    else { s += 2*poly.func(xi); }
  }
  s *= (b-a)/(2*N);
  return s;
}

Esto funciona bien, pero la biblioteca resultante no es muy independiente, ya que necesita que se defina el Modelo de clase en algún lugar. Además, idealmente, el modelo debería poder cambiar de un usuario a otro para que no quisiera arreglarlo en un archivo de encabezado. He intentado usar plantillas de funciones y functores para que esto funcione, pero no es muy independiente, ya que, una vez más, la plantilla debe definirse en un archivo de encabezado (a menos que desee crear una instancia explícita, lo cual no hago).

Entonces, para resumir: ¿hay alguna manera de que mis funciones de integración acepten funciones 1D arbitrarias con un número variable de parámetros de entrada y que sigan siendo lo suficientemente independientes como para compilarlas en una biblioteca independiente? Gracias de antemano por las sugerencias.