double trap(double func(double), double b, double a, double N) {
double j;
double s;
double h = (b-a)/(N-1.0); //Width of trapezia
double func1 = func(a);
double func2;
for (s=0,j=a;j<b;j+=h){
func2 = func(j+h);
s = s + 0.5*(func1+func2)*h;
func1 = func2;
}
return s;
}
The above is my C++ code for a 1D numerical integration (using the extended trapezium rule) of func()
between limits $[a,b]$ using $N-1$ trapezia.
I am actually doing a 3D integration, where this code is called recursively. I work with $N = 50$ giving me decent results.
Other than reducing $N$ further, is anybody able to suggest how to optimise the code above so that it runs faster? Or, even, can suggest a faster integration method?
trapezoidal_integration
instead oftrap
,sum
orrunning_total
instead ofs
(and also use+=
instead ofs = s +
),trapezoid_width
ordx
instead ofh
(or not, depending on your preferred notation for the trapezoidal rule), and changefunc1
andfunc2
to reflect the fact that they are values, not functions. E.g.func1
->previous_value
andfunc2
->current_value
, or something like that. $\endgroup$ – David Z Mar 25 '14 at 1:50