Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Lecture on Numerical Analysis
Dr.-Ing. Michael Dumbser
M. Dumbser
1 / 18
24 / 09 / 2008
Analisi Numerica
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Numerical Integration (Quadrature) of Functions - Motivation
b
Task: compute approximately
f
 f ( x)dx  ?
a
Solution strategy:
• Divide interval [a;b] into n smaller
subintervals
• Approximate f(x) by interpolation
polynomials on the subintervals,
e.g. using Lagrange interpolation
• Integrate these polynomials exactly
on each subinterval and sum up
• So-called Newton-Cotes formulae
h
a
b
n
x
b
xi
n
1
 f ( x)dx    f ( x)dx   h  f ( x  h )d
a
M. Dumbser
2 / 18
i 1 xi 1
i 1
i
0
ba
h
n
xi  a  i  h
i  0,1,...n
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Transformation of the Integration Interval
The computation of numerical quadrature formulae for each sub-interval can be
technically considerably simplified using the following variable substitution:
x  xi  h
xi 1

h  xi 1  xi
dx
h
d
1
1
0
0
~
f ( x )dx  h  f ( xi  h )d : h  f ( )d
xi
Therefore, it is sufficient, without the loss of generality, to consider from now on the
case of numerical integration in the reference interval [0;1].
1

0
M. Dumbser
3 / 18
~
f ( )d  ...
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Newton-Cotes Formulae
The integral of f(x) is approximated in the following steps:
1.
First, the function f(x) is interpolated by a polynomial of degree k inside
each sub-interval. The interpolation points are distributed in an equidistant
manner in each sub-interval.
2.
Second, the interpolation polynomial is integrated analytically.
3.
Steps 1 and 2 produce an approximation of the integral of f(x) in terms
of the function values fi at the interpolation points and the step size h.
M. Dumbser
4 / 18
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Trapezium or Trapezoid Rule
Use linear interpolation polynomials,
i.e. polynomials of degree one in the
subintervals [xi;xi+1]  [0;1]:
f
f    P1    f i
 1
0 1
 f i 1
P1     fi   1  fi 1
1
 P1  d 
b
a
ba
h
n
h
xi  a  i  1h
0
xi 1
x
 P1 x dx 
xi
b
n xi 1
 f ( x)dx  
a
M. Dumbser
5 / 18
i 1
n
1
1
fi  fi 1
2
2
h
 fi  fi 1 
2
h
 f ( xi )  f ( xi 1 )
x P1 ( x)dx  
i 1 2
i
 0
1 0
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Simpson Rule
(Originally Discovered by Johannes Kepler in 1615)
Use quadratic interpolation polynomials,
i.e. polynomials of degree two:
f

  12   1
f    P2    f i

1
0  2 0  1


  0  1
  0  12 
 fi 1
f
 2  0 12  1 i 1 1  01  12 
1
2
b
a
h
ba
n
1
 P2  d 
0
b
h
xi  a  i  1h
6 / 18
a
 f i  1 4 2     f i 1 2 2   
x
2
1
2
1
fi  fi  1  fi 1
6
3 2 6
n xi 1
 f ( x)dx  
M. Dumbser
P2    f i 2 2  3  1 
i 1
xi 1
 P2 x dx 
xi
n


h
f i  4 f i  1  f i 1
2
6

h
f ( xi )  4 f ( xi  1 )  f ( xi 1 )
x P2 ( x)dx  
2
6
i

1
i

Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The 3/8 Rule
Use cubic interpolation polynomials, i.e. polynomials of degree three:


  13   23   1
  0  23   1
f    P3    f i
f

0  13 0  23 0  1 i   13  0 13  23  13  1


  0  13   1
  0  13   23 
 fi 2
 f i 1
2
1 2
 3  0 3  3  3  1
1  01  13 1  23 
1
3
2
3
P3     92 f i  3  2 2  119   92   272 f i  1  3  53  2  23   
3
 272 f i  2  3  43  2  13    92 f i 1  3   2  23  
3
1
 P3  d 
0
b
1
3
3
1
fi  fi  1  fi  2  fi 1
8
8 3 8 3 8
n xi 1
 f ( x)dx  
M. Dumbsera
7 / 18
i 1
n

h
f ( xi )  3 f ( xi  1 )  3 f ( xi  2 )  f ( xi 1 )
x P3 ( x)dx  
3
3
8
i

1
i

Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Isaac Newton
Sir Isaac Newton - (* 4. January 1643 in Woolsthorpe;
† 31. March 1727 in Kensington)
• Physicist, Mathematician, Astronomer and Philosopher
• Together with Leibniz, Newton is one of the inventors of
infinitesimal calculus (differentiation and integration)
• 1667: Fellow of Trinity College, Cambridge
• 1687: Philosophiae Naturalis Principia Mathematica.
Newton discovered the universal law of gravitation
and the laws of motion of classical mechanics
• 1704: Opticks. A corpuscular theory of light.
• From 1703 president of the Royal Society
• Buried in Westminster abbey in London
M. Dumbser
8 / 18
Analisi Numerica
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Gaussian Quadrature Formulae
The previously derived formulae were all of the form
 f  d  w f  
1
0
M
j 1
j
j
j 1
j 
M 1
and were very easy to obtain since an equidistant spacing of the nodes was imposed and
only the weights wj had to be computed. The aim of the Gaussian quadrature formulae is
now to obtain an optimal quadrature formula with a given number of points by making
also the nodes an unknown in the derivation procedure of the quadrature formula and to
come up with an optimal set of nodes xj and weights wj.
An explicit construction strategy for Gaussian integration formulae:
(1) Using M quadrature points, we have M unknowns for the positions and
also M unknowns for the weights, i.e. a total of 2M unknowns.
(2) We need 2M equations to determine uniquely the 2M unknowns.
(3) The equations are obtained by requiring that the integration formula is
exact for polynomials from degree 0 up to degree 2M-1 !
M. Dumbser
9 / 18
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Gaussian Quadrature Formulae
This means we have 2M equations of the form to solve for wj and j.
 P  d  w P  
1
M
i
j 1
0
j i
j  0,1...2 M  1
j
For Pi(), any polynomial of degree i can be used, in particular also the monomial i.
Example 1: One integration point, i.e. M = 1, leading to the two equations:
 P  d   1d  1  w P    w 1
1
1
0
0
0
M
j 1
j 0
j
1
1 M
w j P1  j   w1  1
0 P1  d  0 d  2  
j 1
1
M. Dumbser
10 / 18
1
w1  1, 1 
1
2
1  w1
1
 w11
2
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Gaussian Quadrature Formulae
Example 2: Two integration points, i.e. M = 2, leading to the 4 equations:
1
 1d  1  w  w
1
2
0
1
1
0 d  2  w1  1  w2  2
1
1
3
2
2
2

d



w



w


1
1
2
2

0
1
3

 d 
0
M. Dumbser
11 / 18
1
 w1  13  w2  23
4
1
2
w1  w2  , 1 
1 1
1 1

3,  2  
3
2 6
2 6
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Gaussian Quadrature Formulae
A more efficient and more general way of obtaining the Gaussian quadrature formulae
makes use of so-called orthogonal polynomials Li(), which are the so-called Legendre
polynomials.
First, we define the scalar-product of two functions f and g as follows:
1
f , g   f ( )  g ( )d
0
With this scalar product available, we can define the L2 norm of a function f as
f 
f, f
The set of polynomials Li() is called orthogonal, if it satisfies the relation
M. Dumbser
12 / 18
0
Li , L j  
 i
if
if
i j
i j
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Gaussian Quadrature Formulae
The polynomials Li() can be constructed via Gram-Schmidt orthogonalization
from the monomials M0 = 1, M1 = 2, M3 = 3, … Mn =  n as follows:
We first use the analogy of the scalar product of two functions with the
scalar product already known for vectors:
 
a , b   ai  bi

a 
 
a, a
i
The Gram-Schmidt orthogonalization then proceeds as follows:
M. Dumbser
13 / 18
M1
M1
M 0 : L0
M1,
L0
L0
L0
L0
L1  M 1  M 1 ,
L0
L0
L0
L0
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Gaussian Quadrature Formulae
Instead of performing orthogonalization of vectors, we now perform the
orthogonalization of functions as follows:
L0 ( ) : 1
M 1 , L0
L1 ( )  M 1 ( )  L0
L0 , L0
M 2 , L1
M 2 , L0
L2 ( )  M 2 ( )  L1
 L0
L1 , L1
L0 , L0
M 3 , L2
M 3 , L1
M 3 , L0
L3 ( )  M 3 ( )  L2
 L1
 L0
L2 , L2
L1 , L1
L0 , L0
M. Dumbser
14 / 18
M1
M1
M 0  L0
M1,
L0
L0
L0
L0
L1  M 1  M 1 ,
L0
L0
L0
L0
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Gaussian Quadrature Formulae
Using the orthogonality property of the Legendre polynomials, we find that
1 if i  0
0 Li ( )d  1, Li  0 else
1
(1)
The Gaussian quadrature formulae are written as
1
n
 f ( )d  w
j 1
0
j
f ( j )
(2)
If we now apply formula (2) to the integrals given in (1), we obtain the following
linear equation system for the weights wj, if we suppose the positions j to be known:
1 if i  0
w j Li ( j )  
0 Li ( )d  
j 1
0 else
1
 
Lijw j  ei
M. Dumbser
15 / 18
n
with
Lij  Li ( j )
(3)
(3‘)
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Gaussian Quadrature Formulae
Theorem:
If the n positions j are given be the n roots of the polynomial Ln ( Ln( j ) = 0 ),
and the weights are given by the solution of system (3), then the Gaussian
quadrature rules are exact for polynomials up to degree 2n-1, i.e.
1
n
 p( )d  w p( )
j 1
0
j
j
Proof:
Suppose p() is an arbitrary polynomial of the space of polynomials
of degree 2n-1, i.e. p( )  P2 n 1
Then we can write the polynomial p() as
q( ), r( )  Pn 1
p( )  Ln ( )q( )  r( )
n 1
q( )    k Lk ( )
k 0
M. Dumbser
16 / 18
n 1
r ( )    k Lk ( )
k 0
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Gaussian Quadrature Formulae
Proof (continued):
We have
1
1
 p( )d   L ( )q( )  r( )d 
Ln , q  1, r  0
n
0
0
We also have
w p( )  w L ( )q( )  r( )
n
j 1
n
n
j
j
j 1
n
j
n
j
j
j
n 1
w p( )  w   L ( )
j 1
j
j
j 1
j
n 1
n
k 0
k
k
j
n
w p( )    w L ( )  
j 1
j
j
k 0
k
j 1
This finishes the proof. QED
M. Dumbser
17 / 18
j
k
j
0
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Integration of Improper Integrals
If the integration interval goes to infinity, it can be very useful to change the integration
variables use the following substitution:
dx
1
 2
dt
t
1
1
1
b
a
a
1
1
1
1
 
 
f ( x )dx    f   2 dt   f   2 dt   g (t )dt,
tt
tt
1
1
1
a
b
b
x
b

a
1
1
 t
t
x

ab  0
1 1
g (t )  f    2
t t
Example:
1
0
0
1 1
 f ( x)dx  1 f  t  t 2 dt  1 g t dt
If the integrand is singular at a known position c, than it is usually useful to split the
integral as:
b
c 
b
a
a
c
 f ( x)dx   f ( x)dx   f ( x)dx
M. Dumbser
18 / 18
Note: Gaussian quadrature formulae never use the interval endpoints,
which makes them very useful for the computation of improper integrals!