ME 502 Engineering Computing Lab Page 1 of 10

DEPARTMENT OF MECHANICAL ENGINEERING

INDIAN INSTITUTE OF TECHNOLOGY - GUWAHATI

ME 502 : Engineering Computing Lab

2024-2025, Ist semester

Handout for Lab Session on Numerical

Integration

During numerical simulation, for example in finite element method (FEM),

the need often arises for evaluating the definite integral of a function that has
no explicit antiderivative or whose antiderivative is not easy to obtain. Then,
we need to seek numerical integration techniques in the following situations:

1. Functions do not possess closed form solutions. For example

Z x
2
f (x) = e−t dt (why?) (1)
0

2. Close form solutions exist but these solutions are complex and difficult
to use for calculations.

3. Data for variables are available in the form of a table, but no mathe-
matical relationship between them is known, as is often the case with
experimental data.
Rx
We know that a definite integral of the form for example in 1D x01 f (x) dx
can be treated as the area under the curve y = f (x), enclosed between the
limits x = x0 to x = x1 . If you were to perform the integration manually, one
approach is to superimpose a grid on a graph of the function to be integrated,
and simply count the squares, counting only those covered by 50% or more of
the function, see Figure 1. Provided the grid is sufficiently fine, a reasonably
accurate estimate may be obtained. Figure 1 demonstrates how this may be
achieved.
Although the grid methods and other such graphical approaches can pro-
vide us with rough estimates, they are cumbersome and time-consuming and
the final results are far from satisfactory limits. Also, they are too cumber-
some to be applicable when we go for higher dimensions. A better alternative
approach could be to use a technique that uses simple arithmetic operations to
compute the area. Such an approach can be easily implemented on a computer
to carry out the integration (hence the term numerical integration or numeri-
cal quadrature). Numerical integration techniques are similar to the graphical

Wednesday, 11.09.2024 2024-25, Ist Sem.

 ME 502 Engineering Computing Lab Page 2 of 10

Figure 1: Graphical representation of integral of a function by manual method

by superimposing a grid. The boxes indicated in grey are counted. (source of
the image http:www.damtp.cam.ac.uk/lab/people/sd/lectures/nummeth98/
figure11.htm#FF 11)

or manual methods as both of them use the concept of ”summation” to find

the integral.

The basic method involved in approximating the definite integral - in 1D

Z x2
I(x) = f (x) dx
x1

or in 2D Z x2 Z y2
I(x, y) = f (x, y) dxdy,
x1 y1

or in 3D Z x2 Z y2 Z z2
I(x, y, z) = f (x, y, z) dxdydz
x1 y1 z1

is called numerical integration/numerical quadrature. It uses a sum -

ni
X
wi f (xi )
i=0

in 1D,
nj
ni X
X
wi wj f (xi , yj )
i=0 j=0

in 2D,
nj nk
ni X
X X
wi wj wk f (xi , yj , zk )
i=0 j=0 k=0

Wednesday, 11.09.2024 2024-25, Ist Sem.

ME 502 Engineering Computing Lab Page 3 of 10

in 3D - to approximate Z x2
I(x) = f (x) dx
x1

in 1D, Z x2 Z y2
I(x, y) = f (x, y) dx dy
x1 y1

in 2D, and Z x2 Z y2 Z z2
I(x, y, z) = f (x, y, z) dx dy dz
x1 y1 z1

in 3D. The values xi , yj , or zk are called the sampling points or integration

nodes and the constants wi , wj , and wk are called the weighting coefficients or
simply weights.
There are various methods for selecting the location and number of sam-
pling points. There is a set of methods known as Newton-Cotes rules in which
the sampling points are equally spaced. Another set of methods called Gauss-
Legendre rules, uses sampling points that are not equally spaced, but are de-
signed to provide improved accuracy. In the current exercise we will discuss
these two methods in detail and write computer code to find the solution of
functions numerically. We will also discuss a method known as Romberg inte-
gration that is designed to improve the estimates of Newton-Cotes formulae.
In general, numerical integration methods yield much better results compared
to the numerical differentiation methods (not discussed here ) due to the fact
that the errors introduced in separate subintervals tend to cancel each other.

(A) Newton-Cotes Methods. Newton-Cotes formula is the most popular

and widely numerical integration formula. It forms the basis for a number
of numerical integration methods known as Newton-Cotes methods. The
derivation of Newton-Cotes formula is based on polynomial interpolation.
As pointed out earlier, an nth degree polynomial pn (x) that interpolates
the values of f (x) at n + 1 evenly spaced points can be used to replace
the integrand f (x) of the integral
Z x1
I = f (x) dx (2)
x0

and the resultant formula is called (n + 1) point Newton-Cotes formula.

If the limits of integration x0 and x1 are in the set of interpolating points
xi , i = 0, 1, 2, 3, . . . n, then the formula is referred to as closed form else
open form. Open form formulas are not used for definite integration. In
the present course we only consider the closed form methods. All these
rules can be formulated using either Newton or Lagrange interpolation
polynomial for approximating the function f (x). Some of the Newton-
Cotes methods include

• Trapezoidal Rule: It is the first and the simplest of the Newton-

Cotes methods. Since it is a two-point interpolation formula, it uses
the first order interpolation polynomial p1 (x) for approximating the

Wednesday, 11.09.2024 2024-25, Ist Sem.

 ME 502 Engineering Computing Lab Page 4 of 10

function f (x). This is illustrated in Figure 2. The integral given

by Eq. (2) can be approximated as

Figure 2: Graphical representation of the trapezoidal rule. (source of the

image http:www.damtp.cam.ac.uk/lab/people/sd/lectures/nummeth98/ fig-
ure11.htm#FF 11)

∆x
I ≈ It = (f (x0 ) + f (x1 )). (3)
2
The truncation error in the trapezoidal rule involve f ′′ , so that the
rule gives the exact result when applied to any function whose sec-
ond derivative is identically zero, that is, any polynomial of degree
one or less.
If the range to be integrated is large, the trapezoidal rule can be
improved by dividing the interval x0 , x1 into a number of small
intervals and applying the rule discussed above to each of these
subintervals. The sum of areas of all the subintervals is the integral
of the interval x0 , x1 . This is known as composite trapezoidal rule.
This is illustrated in Figure 3. The integral given by Eq. (2) can
be approximated as
" n−1
#
∆x X
I ≈ Ict = f (x0 ) + 2 f (xi ) + f (x1 ) . (4)
2 i=1

where xi = x0 + i∆x.
• Midpoint Rule: In the midpoint rule the function is evaluated at
1
midpoint of the interval i.e. at xmid = (x0 + x1 ). The integral
2
given by Eq. (2) can be approximated as

I ≈ Imid = ∆xf (xmid ). (5)

This is illustrated in Figure 4 which provides a graphical interpre-

Wednesday, 11.09.2024 2024-25, Ist Sem.

ME 502 Engineering Computing Lab Page 5 of 10

Figure 3: Graphical representation of the com-

pound trapezoidal rule. (source of the image
http:www.damtp.cam.ac.uk/lab/people/sd/lectures/nummeth98/ fig-
ure11.htm#FF 11)

tation of this approach. By evaluating the function f(x) at the

Figure 4: Graphical representation of the midpoint rule. (source of the

image http:www.damtp.cam.ac.uk/lab/people/sd/lectures/nummeth98/ fig-
ure11.htm#FF 11)

midpoint of each interval the error may be slightly reduced relative

to the trapezoidal rule but the method remains of the same order.
Again, as in trapezoidal rule, if the range to be integrated is large,
the accuracy can be increased by dividing the interval x0 , x1 into a
number of small intervals and applying the rule discussed above to
each of these subintervals. This leads to what is called composite
midpoint rule. The integral given by Eq. (2) can be approximated

Wednesday, 11.09.2024 2024-25, Ist Sem.

 ME 502 Engineering Computing Lab Page 6 of 10

as
n−1
X
I ≈ Ic,mid = ∆x f (xi ). (6)
i=0
 
1
where xi = x0 + i + ∆x. This is illustrated in Figure 5 which
2
provides the graphical interpretation of the method. There are two
further advantages of the Mid-point Rule over the trapezoidal rule.
The first is that is requires one fewer function evaluations for a given
number of subintervals, and the second that it can be used more ef-
fectively for determining the integral near an integrable singularity.
The reasons for this are clear from Figure 6.

Figure 5: Composite midpoint rule. (source of the image

http:www.damtp.cam.ac.uk/lab/people/sd/lectures/nummeth98/ fig-
ure11.htm#FF 11)

• Simpson’s 1/3 Rule: Another popular method is Simpsons’s 1/3

rule. Here, the function is approximated by a second-order polyno-
mial which passes through three sampling points. The three points
1
are x = x0 , x = x1 , and x = (x0 + x1 ). Then, the integral
2
given by Eq. (2) can be approximated as

∆x
I ≈ Is, 1 = (f (x0 ) + 4f (xmid ) + f (x1 )) , (7)
3 3
1
where ∆x = (x1 − x0 ). It is important to note that the Simp-
2
sons’s 1/3 rule is exact up to degree 3, although it is based on
quadratic equation. Similar to the composite trapezoidal rule,we
can construct a composite Simpson’s 1/3 rule to improve the es-
timate of the integral. Here again, the integration interval is di-
vided into n number of segments of equal width, where n is an
1
even number. Them the step size is ∆x = (x1 − x0 ). Again
n
Wednesday, 11.09.2024 2024-25, Ist Sem.
ME 502 Engineering Computing Lab Page 7 of 10

Figure 6: Applying the midpoint rule where the singular inte-

grand would cause the trapezoidal rule to fail. (source of the
image http:www.damtp.cam.ac.uk/lab/people/sd/lectures/nummeth98/ fig-
ure11.htm#FF 11)

xi = x0 + i ∆x, where i = 0, 1, 2, . . . n. Then, the integral given

by Eq. (2) can be approximated as
 
n/2 (n/2)−1
∆x  X X
I ≈ Isc, 1 = f (x0 ) + 4 f (x2i−1 ) + 2 f (x2i ) + f (x1 ) ,
3 3 i=1 i=1
(8)
• Simpson’s 3/8 Rule: Simpson’s 1/3 rule was derived using three
sampling points that fit a quadratic equation. We can extend this
approach to incorporate four sampling points so that the rule can
be exact for f (x) of degree 3. Then, the integral given by Eq. (2)
can be approximated as
3∆x
I ≈ Is, 3 = (f (x0 ) + 3f (xm ) + 3f (xn ) + f (x1 )) , (9)
8 8
1
where ∆x = (x1 − x0 ). This equation is known as the Simpson’s
3
3/8 rule. Here, xm = x0 + ∆x and xn = x0 + 2 ∆x.
• Boole’s Rule: This method is based on five sampling points. Then,
the integral given by Eq. (2) can be approximated as

3∆x
I ≈ Ib = (7f (x0 ) + 32f (xm ) + 12f (xn ) + 32f (xo ) + 7f (x1 )) ,
45
(10)
1
where ∆x = (x1 − x0 ). Here, xm = x0 + ∆x, xn = x0 + 2 ∆x,
4
and xo = x0 + 3 ∆x.

(B) Gauss-Legendre Integration Methods or Gauss quadrature rule.

Wednesday, 11.09.2024 2024-25, Ist Sem.

 ME 502 Engineering Computing Lab Page 8 of 10

We have discu–ssed so far a set of rules based on Newton-Cotes formula.

In Newton-Cotes formula all the sampling points were evenly placed
within the range of the integral. Gauss-Legendre methods (or popularly
known simply as Gauss integration method or Gauss quadrature rule) is
based on the concept that the accuracy of the numerical integration can
be improved by choosing the sampling points wisely, rather than on the
basis of equal spacing. In this method the sampling points lie inside the
domain of integration and the end points are not used in the computa-
tion. Then, the integral given by Eq. (2) and written for x0 = −1 and
x1 = 1, x = ξ can be approximated as
Z +1 n = nξ
X
I = f (ξ)dξ ≈ Ig = wi f (ξi ), (11)
−1 i=1

where ξi are the sampling points (better known as Gauss points) and
wi = w(ξi ) is called the Gauss point weight corresponding to the Gauss
point ξi . Here, n is the called the order of the Gauss integration method.
Hence, Eq. (11) is called the n-point formula. The problem is to com-
pute the locations of the sampling points and to choose the appropriate
”weights” wi for each Gauss point. Gauss integration method assumes
that x0 = −1 and x1 = +1. In doing so the Gauss point location and
weights are be readily available1 . A MatLab script file lgwt.m is available
on the web which produces the Legendre-Gauss weights and nodes for
computing the definite integral of a continuous function on some interval
(x0 , x1 ) 2 . Table below shows the Gauss point coordinates (locations)
and weight for up to n = 3.

n i or j or k wi or wj or wk ξi or ηj or ζk
1 1 2 0
2 1 1.00000 -0.57735
2 1.00000 +0.57735
3 1 0.55556 +0.77460
2 0.88889 0.00000
3 0.55556 -0.77460

Considering polynomial integrands, it has been proven that the use of n

Gauss points gives the exact result of a polynomial integrand of up to the
order of p = 2n − 1. For example, if the integrand is a linear function,
we have 2n − 1 = 1, which gives n = 1. This means that for a linear
integrand p = 1, one Gauss point will be sufficient to give the exact
result. of the integration. If the integrand is a polynomial of order p = 3,
then we have 2n − 1 = 3, which gives n = 2. This means that for an
integrand of a third order polynomial, the use of two Gauss points will be
1
see https://pomax.github.io/bezierinfo/legendre-gauss.html for values upto n = 64.
2
https://in.mathworks.com/matlabcentral/fileexchange/4540-legendre-gauss-
quadrature-weights-and-nodes

Wednesday, 11.09.2024 2024-25, Ist Sem.

ME 502 Engineering Computing Lab Page 9 of 10

sufficient to give the exact result. The use of more than two Gauss points
will still give the same results, but takes more computational time. The
Gauss integration procedure is by no means limited to one dimensions.
In finite element analysis, for example, integrals of the form
Z +1 Z +1
I = f (ξ, η)dξdη (12)
−1 −1
Z +1 Z +1 Z +1
I = f (ξ, η, ζ)dξdηdζ (13)
−1 −1 −1

are frequently encountered. For such cases the Gauss integration is sam-
pled in two or three directions, as follows
Z +1 Z +1 nξ nη
X X
f (ξ, η)dξdη ≈ wξ wη f (ξi , ηj ) (14)
−1 −1 i=1 j=1
Z +1 Z +1 nξ nη nζ
XXX
f (ξ, η, ζ)dξdηdζ ≈ wξ w(ξi ) wη wζ f (ξi , ηj , ζ(15)
k)
−1 −1 i=1 j=1 k=1

Here, nξ , nη , and nζ are the number of Gauss points in the ξ, η, and ζ

directions. Also, wξ = w(ξi ), wη = w(ηj ) and wζ = w(ζj ) are Gauss
points weights corresponding to Gauss points ξi , ηj , and ζk respectively.
The values can be read from the table above once the number of Gauss
points in each direction is fixed.
Suggested References
1. M. T. Heath, Scientific Computing - An Introductory Survey, Revised Sec-
ond Edition, SIAM, 2018 (also see: https://heath.cs.illinois.edu/scicomp/).
2. R. L. Burden and J. D. Faires, Numerical Analysis, Ninth Edn., Thomson
Brooks/Cole, 2011.
3. K. E. Atkinson, An Introduction to Numerical Analysis, Second Edn., Wi-
ley, 2008.
4. S. D. Conte and C. de Boor, Elementary Numerical Analysis, Third Edition,
Tata McGraw-Hill Education, 2005.
5. J. D. Hoffman, Numerical Methods for Engineers and Scientists, Second
Edn., CRC Press, 2001.
6. G. M. Phillips and P. J. Taylor, Theory and Applications of Numerical
Analysis, Second Edition, Academic Press, 1996.
7. F.B. Hildebrand, Introduction to Numerical Analysis, Second (Revised)
Edition, Courier Dover Publications, 1987.
8. R. W. Hamming, Numerical Methods for Scientists and Engineers, Second
Edition, Dover, 1986

Wednesday, 11.09.2024 2024-25, Ist Sem.

 ME 502 Engineering Computing Lab Page 10 of 10

IMPORTANT INSTRUCTIONS
(a) Please write a separate function in C for each integration rule. Also, make
a dedicated function for polynomial evaluation at the sampling point.
This will help you in using the same programme for different polynomials
by just rewriting the function only without the need to change the whole
code.

(b) During the evaluation of this assignment, for each problem please report
the error with respect to the exact answer, where ever applicable. So, it
is preferred that you make tabulate all your results in a single sheet.

(c) Also, it will be desirable that you also carry out a parametric study like
increasing the number of intervals, or increasing the number of Gauss
integration points and report your observation.
Rx
Question 1. Evaluate the integral I = x01 (x3 + 1)dx using the trapezoidal
rule for the intervals (a) (1, 2) and (b) (1, 1.5). Also compare the results
with the exact answer in each case. (Note: Make a function by the name
trapezoidal in C for this.)
R1
Question 2. Compute the integral I = −1 ex dx using the composite trape-
zoidal rule for (a) n = 2 and (b) n = 4. Compare your result with
exact result Iexact = 2.35040. (Note: Make a function by the name
trapezoidalC in C for this.)
R 10
Question 3. Compute the integral I = 0 1 − e−x/2 dx using the com-

posite trapezoidal rule for ∆x = 0.5. Compare your result with exact
result. (Note: Make a function by the name trapezoidalC in C for this.)

Question 4. Evaluate the integral in Question 1R using using Simpson’s 1/3

π√
rule and Simpson 3/8 rule while the integral 0 sinx dx using composite
Simpson’s 1/3 rule. (Note: Make a function by the name simpson13,
simpson38, simpson13C in C for this.)
R +1
Question 5. Evaluate the integral −1 (ξ 2 − 3ξ + 7) dξ using Gauss quadra-
ture rule so that the result is exact. (Note: Make a function by the name
gauss1D in C for this.)
R +1 R +1
Question 6. Evaluate the integral −1 −1 (ξ 3 − 1) (η − 1)2 dξ dη using Gauss
quadrature rule so that the result is exact. (Note: Make a function by
the name gauss2D in C for this.)
R +1 R +1 R +1
Question 7. Evaluate the integral −1 −1 −1 ξ 2 (η 2 − 1) (ζ 4 − 2) dξ dη dζ
using Gauss quadrature rule so that the result is exact. (Note: Make a
function by the name gauss3D in C for this.)

NOTE: In all the above examples make a function function evaluation

for computing the function at the sampling points.

Wednesday, 11.09.2024 2024-25, Ist Sem.