Computational Physics - Department of Physics

5.1 Newton-Cotes Quadrature 113

in the calling function. We note thata,bandnare called by value, whiletrapez_sumand
the user defined functionMyFunctionare called by reference.
The name trapezoidal rule follows from the simple fact that it has a simple geometrical
interpretation, it corresponds namely to summing up a series of trapezoids, which are the
approximations to the area below the curvef(x).
Another very simple approach is the so-called midpoint or rectangle method. In this case
the integration area is split in a given number of rectangleswith lengthhand height given
by the mid-point value of the function. This gives the following simple rule for approximating
an integral

I=

∫b

a

f(x)dx≈h

N

∑

i= 1

f(xi− 1 / 2 ), (5.5)

wheref(xi− 1 / 2 )is the midpoint value offfor a given rectangle. We will discuss its truncation
error below. It is easy to implement this algorithm, as shownhere

http://folk.uio.no/mhjensen/compphys/programs/chapter05/cpp/rectangle.cpp
doublerectangle_rule(doublea,doubleb,intn,double(func)(double))
{
doublerectangle_sum;
doublefa, fb, x, step;
int j;
step=(b-a)/((double) n);
rectangle_sum=0.;
for(j = 0; j <= n; j++){
x = (j+0.5)step+;// midpoint of a given rectangle
rectangle_sum+=(func)(x);// add value of function.
}
rectangle_sum= step;// multiply with step length.
returnrectangle_sum;
}// end rectangle_rule

The correct mathematical expression for the local error forthe rectangular ruleRi(h)for
elementiis ∫
h
−h

f(x)dx−Ri(h) =−h

3

24

f(^2 )(ξ),

and the global error reads
∫b
a

f(x)dx−Rh(f) =−b−a

24

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

whereRhis the result obtained with rectangular rule andξ∈[a,b].
Instead of using the above first-order polynomials approximations forf, we attempt at
using a second-order polynomials. In this case we need threepoints in order to define a
second-order polynomial approximation

f(x)≈P 2 (x) =a 0 +a 1 x+a 2 x^2.

Using again Lagrange’s interpolation formula we have

P 2 (x) =

(x−x 0 )(x−x 1 )

(x 2 −x 0 )(x 2 −x 1 )

y 2 +

(x−x 0 )(x−x 2 )

(x 1 −x 0 )(x 1 −x 2 )

y 1 +

(x−x 1 )(x−x 2 )

(x 0 −x 1 )(x 0 −x 2 )

y 0.

Inserting this formula in the integral of Eq. (5.2) we obtain
∫+h
−h

f(x)dx=h

3

(fh+ 4 f 0 +f−h)+O(h^5 ),