Computational Physics - Department of Physics

112 5 Numerical Integration

Hereafter we use the shorthand notationsf−h=f(x 0 −h),f 0 =f(x 0 )andfh=f(x 0 +h). The
correct mathematical expression for the local error for thetrapezoidal rule is

∫b

a

f(x)dx−

b−a

2 [f(a)+f(b)] =−

h^3

12 f

( 2 )(ξ),

and the global error reads
∫b
a

f(x)dx−Th(f) =−

b−a

12

h^2 f(^2 )(ξ),

whereThis the trapezoidal result andξ∈[a,b].
The trapezoidal rule is easy to implement numerically through the following simple algo-
rithm

• Choose the number of mesh points and fix the step.


• calculatef(a)andf(b)and multiply withh/ 2


• Perform a loop overn= 1 ton− 1 (f(a)andf(b)are known) and sum up the terms
  f(a+h) +f(a+ 2 h) +f(a+ 3 h) +···+f(b−h). Each step in the loop corresponds to a
  given valuea+nh.


• Multiply the final result byhand addh f(a)/ 2 andh f(b)/ 2.

A simple function which implements this algorithm is as follows

http://folk.uio.no/mhjensen/compphys/programs/chapter05/cpp/trapezoidal.cpp
doubletrapezoidal_rule(doublea,doubleb,intn,double(func)(double))
{
doubletrapez_sum;
doublefa, fb, x, step;
int j;
step=(b-a)/((double) n);
fa=(func)(a)/2. ;
fb=(func)(b)/2. ;
TrapezSum=0.;
for(j=1; j <= n-1; j++){
x=jstep+a;
trapez_sum+=(func)(x);
}
trapez_sum=(trapez_um+fb+fa)step;
returntrapez_sum;
}// end trapezoidal_rule

The function returns a new value for the specific integral through the variabletrapez_sum.
There is one new feature to note here, namely the transfer of auser defined function called
funcin the definition

voidtrapezoidal_rule(doublea,doubleb,intn,double*trapez_sum,

double(*func)(double) )

What happens here is that we are transferring a pointer to thename of a user defined func-
tion, which has as input a double precision variable and returns a double precision number.
The functiontrapezoidal_ruleis called as

trapezoidal_rule(a, b, n, &MyFunction )