long time browser, first time asker here. I've written a number of scripts for doing various 1D numerical integration methods and compiled them into a library. I would like that library to be as flexible as possible regarding what it is capable of integrating.
Here I include an example: a very simple trapezoidal rule example where I pass a pointer to the function to be integrated.
// 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;
}
This works great for simple functions that only take one argument. Example:
double a = trap(sin,0,1);
However, sometimes I may want to integrate something that has more parameters, like a quadratic polynomial. In this example, the coefficients would be defined by the user before the integration. Example code:
// arbitrary quadratic polynomial
double quad(double A, double B, double C, double x) {
return (A*pow(x,2) + B*x + C);
}
Ideally, I would be able to do something like this to integrate it:
double b = trap(quad(1,2,3),0,1);
But clearly that doesn't work. I have gotten around this problem by defining a class that has the coefficients as members and the function of interest as a member function:
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); }
};
However, then my integration function needs to change to take an object as input instead of a function pointer:
// 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;
}
This works fine, but the resulting library is not very independent, since it needs the class Model to be defined somewhere. Also, ideally the Model should be able to change from user-to-user so I wouldn't want to fix it in a header file. I have tried to use function templates and functors to get this to work but it is not very independent since again, the template should be defined in a header file (unless you want to explicitly instantiate, which I don't).
So, to sum up: is there any way I can get my integration functions to accept arbitrary 1D functions with a variable number of input parameters while still remaining independent enough that they can be compiled into a stand-alone library? Thanks in advance for the suggestions.
What you need is templates and std::bind()
(or its boost::bind()
counterpart if you can't afford C++11). For instance, this is what your trap()
function would become:
template<typename F>
double trap(F&& f, 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;
}
Notice, that we are generalizing from function pointers and allow any type of callable objects (including a C++11 lambda, for instance) to be passed in. Therefore, the syntax for invoking the user-provided function is not *f(param)
(which only works for function pointers), but just f(param)
.
Concerning the flexibility, let's consider two hardcoded functions (and pretend them to be meaningful):
double foo(double x)
{
return x * 2;
}
double bar(double x, double y, double z, double t)
{
return x + y * (z - t);
}
You can now provide both the first function directly in input to trap()
, or the result of binding the last three arguments of the second function to some particular value (you have free choice on which arguments to bind):
#include <functional>
int main()
{
trap(foo, 0, 42);
trap(std::bind(bar, std::placeholders::_1, 42, 1729, 0), 0, 42);
}
Of course, you can get even more flexibility with lambdas:
#include <functional>
#include <iostream>
int main()
{
trap(foo, 0, 42);
trap(std::bind(bar, std::placeholders::_1, 42, 1729, 0), 0, 42);
int x = 1729; // Or the result of some computation...
int y = 42; // Or some particular state information...
trap([&] (double d) -> double
{
x += 42 * d; // Or some meaningful computation...
y = 1; // Or some meaningful operation...
return x;
}, 0, 42);
std::cout << y; // Prints 1
}
And you can also pass your own stateful functors tp trap()
, or some callable objects wrapped in an std::function
object (or boost::function
if you can't afford C++11). The choice is pretty wide.
Here is a live example.
What you trying to do is to make this possible
trap( quad, 1, 2, 3, 0, 1 );
With C++11 we have alias template and variadic template
template< typename... Ts >
using custom_function_t = double (*f) ( double, Ts... );
above define a custom_function_t
that take a double and variable numbers of arguments.
so your trap
function becomes
template< typename... Ts >
double trap( custom_function_t<Ts...> f, Ts... args, 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, args...); }
else { s += 2*f(xi, args...); }
}
s *= (b-a)/(2*N);
return s;
}
Usage:
double foo ( double X ) {
return X;
}
double quad( double X, double A, double B, double C ) {
return(A*pow(x,2) + B*x + C);
}
int main() {
double result_foo = trap( foo, 0, 1 );
double result_quad = trap( quad, 1, 2, 3, 0, 1 ); // 1, 2, 3 == A, B, C respectively
}
Tested on Apple LLVM 4.2 compiler.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With