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 / 16
01 / 12 / 2008
Analisi Numerica
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Linear Algebra
A standard task in linear algera is finding the solution of a system of linear equations of the form
A X  B
The solution to this task is given by the so-called Gauß-algorithm:
Define a new matrix C as
C  Cij   A, B 
Now the matrix C is modified by a sequence of operations on its rows to transform its left
part into the unit matrix.
The admissible row-operations are:
(1) addition / subtraction of multiples of rows
(2) multiplication / division of a row by a constant factor
(3) exchange of rows
When the left part of the matrix C has been transformed into the unit matrix, the right part contains
the solution X of the equation system.
M. Dumbser
2 / 16
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Gauß-Algorithm I
In detail, the Gauß-algorithm for the N x N matrix A and the N x M matrix B proceeds as follows:
Part 1 is a forward loop, which transforms the matrix A into a normalized upper triangular form.
Loop over all rows. The loop index i runs from 1 to N
(1) check the coefficient C(i,i) on the diagonal of the row.
(2) IF it is zero, look for the next row j with row index j > i which has a
nonzero element in column i.
IF there is no such row, then the matrix is singular. EXIT.
ELSE exchange rows i and j in the matrix C.
(3) Divide the whole row i by C(i,i). Its first non-zero entry is now equal to 1.
(4) Subtract a suitable multiple of row i from all rows j > i in order to eliminate all entries in
column i of those rows j.
end loop
M. Dumbser
3 / 16
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Gauß-Algorithm II
In detail, the Gauß-algorithm for the N x N matrix A and the N x M matrix B proceeds as follows:
Part 2 is a very simple backward loop, which transforms the matrix A into the unit matrix.
Loop over all rows. The loop index i runs from N to 1
Subtract a suitable multiple of row i from all rows j < i in order to eliminate all
entries in column i of those rows j.
end loop
All operations have been performed on the augmented matrix C. The type of row-operations is
completely determined by the left part if matrix C, i.e. by the original matrix A. The right part, i.e.
matrix B, is transformed with the same operations in a passive way.
At the end of the Gauß-algorithm, the right part of the matrix C, i.e. where B was located
initially, we find the solution X of the equation system.
To compute the inverse of a matrix using the Gauß-algorithm, the right hand side B must
be simply set to the the N x N unit matrix, since we have:
M. Dumbser
4 / 16
A X  I

1
XA
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Exercise 1
Write a MATLAB function gauss.m that accepts a general N x N matrix A and a general
N x M matrix B as input argument and that returns the N x M matrix X
that solves the matrix equation
A X  B
(1) Use the Gauß algorithm to solve the equation system
1 2 0 4 
1


 
 0  2 0  4   1
 0 1 0 1   x  1
2
2


 
1
1
1
1 
1
4
2 
2 8
(2) Use the Gauß algorithm to compute the inverse of a random 5 x 5 matrix A
generated with the MATLAB command rand as follows:
A = rand(5).
Hint: to generate a 5x5 unit matrix, use the following MATLAB command:
B = eye(5).
M. Dumbser
5 / 16
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
A Special Case of Gauß-Eliminiation: The Thomas-Algorithm
In the case where the matrix A is tridiagonal, i.e. with the special structure
 b1

 a2
0
A
0
0

0

c1
0
0
0
b2
c2
0
0
a3 b3
c3
0
0
ai
bi
ci
0
0
0
0
a N 1 bN 1
0
aN
0 

0 
0 

0 
cN 1 

bN 
there is a very fast and efficient special case of the Gauß-algorithm, called the Thomas-algorithm.
Since it is a variation of the Gauß-algorith, the Thomas algorithm is a direct method to solve general
linear tridiagonal systems. As the original Gauß-algorithm, it proceeds in two stages, one forward
elimination and one back-substitution.
M. Dumbser
6 / 16
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
The Thomas-Algorithm
Part I: Forward elimination.
c1 : c1 / b1; d1 : d1 / b1
 : 1 / bi  ci 1ai 
ci : ci  
d i : d i  ai  d i 1   
Part II: Back substitution
xN  d N
xi  d i  ci  xi 1
M. Dumbser
7 / 16
Analisi Numerica
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Exercise 2
Write a MATLAB function Thomas.m that realizes the Thomas algorithm and that
accepts as input four vectors of equal length a, b, c and d, where a contains the
lower diagonal elements, b contains the diagonal elements, c contains the upper
diagonal elements and the vector d contains the right hand side of the linear
equation system Ax = d. The output of Thomas.m is the solution vector x that
satisfies the above-mentioned system.
To test your code, proceed in the following steps:
(1) Generate the vectors a,b,c and d randomly using the MATLAB command rand.
(2) Assemble the full matrix A using the vectors a, b and c.
(3) Solve the system Ax = d directly using MATLAB.
(4) Solve the sytem Ax = d using the function call x = thomas(a,b,c,d).
M. Dumbser
8 / 16
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Conjugate Gradient Method
For large matrices A, the number of operations of the Gauß algorithm grows with N3. For large but
sparse matrices (many entries are zero), it may be more efficient to use an iterative scheme for the
solution of the system. For symmetric, positive definite matrices A, a very efficient algorithm can be
constructed, which is the so-called conjugate gradient method, which we have already seen for
nonlinear, multi-variate optimization.
We start with some definitions:
(1) A matrix is symmetric positive definite, if
A A
T
and
 

x T A x  0 x  0
(2) Two vectors u and v are called conjugate with respect to matrix A if
 
uT A v  0
(3) If we knew N conjugate directions pj (j = 1 .. N), we could write the solution x of




x  1 p1  2 p2  ...   N pN
The coefficients  are given by:
T 
p b
 j  T j 
pj Apj
M. Dumbser
9 / 16
(4) The term
 

r  b  A x is called the residual
 
A x  b
as
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Conjugate Gradient Method
The idea of the CG method is to minimize the functional
 1 T  T 
g ( x )  x Ax  b x
2
Its gradient is:
 

g

A
x

b


r

x
Using the steepest descent method in direction p, we obtain the minimum of g starting from xk as follows:
 

x  xk   p
g


pT
 T
 x  g  T  
      p Ax  b  0
   x

 
AxK  p   b  0





T 

T 
p b  AxK
p rk


 
 
p T Ap
p T Ap
M. Dumbser
10 / 16

Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Conjugate Gradient Method
A preliminary version of the conjugate gradient method now reads as

x0  0
DO k = 0... N - 1
 

r0  b  Ax0
 
p0  r0
From definition (3) we know that the algorithm will terminate
with the exact solution after at most N iterations


vk  Apk
 T
pk rk
k   T 
1D minimization along direction pk
pk v k



Move point x into the minimum along direction pk
xk 1  xk  k pk






rk 1  b  Axk 1  b  Axk  k vk
Compute new residual ( = negative gradient



= steepest descent direction)
rk 1  rk  k vk
 
k
p


j Ark 1 
Do not move in the steepest descent direction,
pk 1  rk 1     p j
but compute a new conjugate direction pk+1 that
j 0 p j A p j
takes into account all the previous search

ENDDO
M. Dumbser
11 / 16

directions. As in nonlinear minimization, this
guarantees faster convergence.
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Conjugate Gradient Method
After some algebra, the final conjugate gradient method is given as follows

x0  0
 

r0  b  Ax0
DO k = 0... N - 1


vk  Apk
k
k   T 
pk v k



xk 1  xk  k pk



rk 1  rk  k vk
 
 k 1  rkT1rk 1


 k 1 
pk 1  rk 1 
pk
k
ENDDO
M. Dumbser
12 / 16
 
p0  r0
 
0  r0T r0
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The GMRES Method of Saad and Schultz (1986)

 
The GMRES method of Saad and Schultz minimizes the residual norm g ( x )  Ax  b
 
x  x0  0
 

r0  b  Ax0

 1  r0


 
v1  r0 /  1 m1  Av1
for j = [1:N]
 
i=1...j
hij  viT m j j
 

w  m j   h (i, j ) vi
 i 1
h j 1 j  w
for
i = [1:j-1]
q  si 1hij  ci 1hi 1 j
hij  p
end
  h 2jj  h 2j 1 j
h jj : 
 j 1 : s j 1 j
if( j+1 > tol )
M. Dumbser
13 / 16

p  ci 1hij  si 1hi 1 j
hi 1 j  q
s j 1  h j 1 j / 
 j : c j 1 j
(see next page)
c j 1  h jj / 
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The GMRES Method of Saad and Schultz (1986)
...contiuation of the previous page
for j=[1:N]

if(j+1 > tol)


v j 1  w / h j 1 j


m j 1  Av j 1
else
for i=[j:-1:1]
i   i 
j
h 
k i 1
 i   i / hii
end
j
 

x  x0    i vi
end
M. Dumbser
14 / 16
end
i 1
ik
k
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Exercise 3
Write a MATLAB function CG.m that solves the system
 
A x  b
The function should check if A verifies the necessary symmetry
property for the CG method.
Solve the system with
2 1 0 


A  1 4 5 
 0 5 10 


M. Dumbser
15 / 16
1
  
b   2
 3
 
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Linear Least-Squares Method
In the case where we have more equations than unknowns in the linear system
 
A x  b
 
A x  b  0
we call the equation system an overdetermined system. Typical applications for such
overdetermined systems can be found in data analysis (linear regression) where so-called
best-fit curves have to be computed e.g. from observations or experimental data.
In this case, the matrix is no longer a N x N square matrix, but a M x N rectangular matrix
with M > N. In general, such systems do not have a solution in the original sense of the
equation. However, according to ideas of Adrien Marie Legendre and Carl Friedrich Gauß,
we can try to find the best possible solution of the system by minimizing the following goal
function:

 


 
 T
  2
g( x)  A  x  b  A  x  b  A  x  b  r
M. Dumbser
16 / 16
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Linear Least-Squares Method
The minimum of g(x) can be computed analytically using the differential criterion for a minimum:
 T T  T  T T  T 
g( x)  x A A x  b A x  x A b  b b


g
T
T
  2A A x  2A b  0
x
The least-squares solution of the overdetermined equation system can therefore be found by solving
the so-called normal equation


T
A Ax  A b
T


T
x A A
A A 
T
M. Dumbser
17 / 16
1
T
A

1

A b
T
is called the pseudo – inverse of the N x M matrix A.
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
The Linear Least-Squares Method
M. Dumbser
18 / 16
Analisi Numerica
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Exercise 4
Write a MATLAB script LSQ.m that computes the coefficients of the parabola for an
oblique throw from position (x0,y0) with horizontal and vertical speed (vx,vy). The
gravitation constant is g = 9.81.
The „measurements“ are available in form of points (xi,yi) that are generated from
the physical parabola of the throw plus a random (positive and negative) perturbation
of amplitude A.
Use the linear Least-Squares Method to solve the problem, solving directly the
normal equation of the problem.
Plot the „measured data“ as well as the parabola obtained from the least squares
method in the same figure.
Compare the exact coefficients of the parabola obtained from the physics of an
oblique throw with the coefficients obtained by the least squares method in function
of the number of points in the measurements.
M. Dumbser
19 / 16
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Carl Friedrich Gauß
German mathematician, astronomer and physicist
30/04/1777 – 23/02/1855
Significant contributions to
• Linear algebra (Gauß algorithm, Least Squares Method)
• Statistics (Gaussian normal probability distribution)
• Potential theory and differential calculus (Gauß
divergence theorem)
• Numerical analysis (Gaussian quadrature rules)
Already during his life, he was called the „prince of mathematics“.
He produced more than 250 scientific contributions to his research fields.
M. Dumbser
20 / 16
Analisi Numerica