Applied Linear Algebra and Optimization Using MATLAB
Applied Linear Algebra and Optimization Using MATLAB
and Optimization
using MATLAB
R
License, Disclaimer of Liability, and Limited Warranty
By purchasing or using this book (the “Work”), you agree that this license grants
permission to use the contents contained herein, but does not give you the right
of ownership to any of the textual content in the book or ownership to any of
the information or products contained in it. This license does not permit up-
loading of the Work onto the Internet or on a network (of any kind) without the
written consent of the Publisher. Duplication or dissemination of any text, code,
simulations, images, etc. contained herein is limited to and subject to licensing
terms for the respective products, and permission must be obtained from the
Publisher or the owner of the content, etc., in order to reproduce or network
any portion of the textual material (in any media) that is contained in the Work.
The author, developers, and the publisher of any accompanying content, and
anyone involved in the composition, production, and manufacturing of this work
will not be liable for damages of any kind arising out of the use of (or the inabil-
ity to use) the algorithms, source code, computer programs, or textual material
contained in this publication. This includes, but is not limited to, loss of revenue
or profit, or other incidental, physical, or consequential damages arising out of
the use of this Work.
The sole remedy in the event of a claim of any kind is expressly limited to
replacement of the book, and only at the discretion of the Publisher. The use of
“implied warranty” and certain “exclusions” vary from state to state, and might
not apply to the purchaser of this product.
Applied Linear Algebra
and Optimization
using MATLAB
R
This publication, portions of it, or any accompanying software may not be reproduced
in any way, stored in a retrieval system of any type, or transmitted by any means,
media, electronic display or mechanical display, including, but not limited to,
photocopy, recording, Internet postings, or scanning, without prior permission in
writing from the publisher.
ISBN: 978-1-9364200-4-9
The publisher recognizes and respects all marks used by companies, manufacturers,
and developers as a means to distinguish their products. All brand names and
product names mentioned in this book are trademarks or service marks of their
respective companies. Any omission or misuse (of any kind) of service marks or
trademarks, etc. is not an attempt to infringe on the property of others.
1112133 2 1
Our titles are available for adoption, license, or bulk purchase by institutions,
corporations, etc. For additional information, please contact the
Customer Service Dept. at 1-800-758-3756 (toll free).
Preface xv
Acknowledgments xix
Appendices 917
Bibliography 1129
Index 1145
Preface
The main approach used in this book is quite different from currently
available books, which are either too theoretical or too computational. The
approach adopted in this book lies between the above two extremities. The
book fully exploits MATLAB’s symbolic, numerical, and graphical capabil-
ities to develop a thorough understanding of linear algebra and optimiza-
tion algorithms.
The book covers two distinct topics: linear algebra and optimization
theory. Linear algebra plays an important role in both applied and
theoretical mathematics, as well as in all of science and engineering,
xvi Preface
ing with the intent of this book, this chapter presents the mathematical
formulations of basic linear programming problems. In Chapter 7, we de-
scribe nonlinear programming formulations. We discuss many numerical
methods for solving unconstrained and constrained problems. In the be-
ginning of the chapter some of the basic mathematical concepts useful in
developing optimization theory are presented. For unconstrained optimiza-
tion problems we discuss the golden-section search method and quadratic
interpolation method, which depend on the initial guesses that bracket the
single optimum, and Newton’s method, which is based on the idea from
calculus that the minimum or maximum can be found by solving f 0 (x) = 0.
For the functions of several variables, we use the steepest descent method
and Newton’s method. For handling nonlinear optimization problems with
constraints, we discuss the generalized reduced-gradient method, Lagrange
multipliers, and KT conditions. At the end of the chapter, we also discuss
quadratic programming problems and the separable programming prob-
lems.
My sincere thanks are also due to the Deanship of the Scientific Re-
search Center, College of Science, King Saud University, Riyadh, KSA,
for financial support and for providing facilities throughout the research
project No. (Math/2008/05/B).
It has taken me five years to write this book and thanks must go to my
long-suffering family for my frequent unsocial behavior over these years. I
am profoundly grateful to my wife Saima, and our children Fatima, Usman,
xx Acknowledgments
1.1 Introduction
When engineering systems are modeled, the mathematical description is
frequently developed in terms of a set of algebraic simultaneous equations.
Sometimes these equations are nonlinear and sometimes linear. In this
chapter, we discuss systems of simultaneous linear equations and describe
the numerical methods for the approximate solutions of such systems. The
solution of a system of simultaneous linear algebraic equations is proba-
bly one of the most important topics in engineering computation. Prob-
lems involving simultaneous linear equations arise in the areas of elasticity,
electric-circuit analysis, heat transfer, vibrations, and so on. Also, the
numerical integration of some types of ordinary and partial differential
equations may be reduced to the solution of such a system of equations. It
has been estimated, for example, that about 75% of all scientific problems
require the solution of a system of linear equations at one stage or another.
It is therefore important to be able to solve linear problems efficiently and
accurately.
1
2 Applied Linear Algebra and Optimization using MATLAB
For example,
4x1 − 2x2 = 5
3x1 + 2x2 = 4
is a system of two equations in two variables x1 and x2 , and
There are three possible types of linear systems that arise in engineering
problems, and they are described as follows:
1. If there are more equations than unknown variables (m > n), then
the system is usually called overdetermined. Typically, an overdeter-
mined system has no solution. For example, the following system
4x1 = 8
3x1 + 9x2 = 13
3x2 = 9
has no solution.
2. If there are more unknown variables than the number of the equations
(n > m), then the system is usually called underdetermined. Typi-
cally, an underdetermined system has an infinite number of solutions.
For example, the system
x1 + 5x2 = 45
3x2 + 4x3 = 21
4 Applied Linear Algebra and Optimization using MATLAB
2x1 + x2 − x3 = 1
x1 − 2x2 + 3x3 = 4
x1 + x2 = 1
x1 = 1, x2 = 0, x3 = 1.
5x1 + x2 + x3 = 4
3x1 − x2 + x3 = −2
x1 + x2 = 3
does not have a unique solution since the equations are not linear
independent; the first equation is equal to the second equation plus
twice the third equation.
Matrices and Linear Systems 5
Every system of linear equations has either no solution, exactly one solu-
tion, or infinitely many solutions. •
x1 + x2 = 1
2x1 + 2x2 = 3.
From the graphs (Figure 1.1(a)) of the given two equations we can see
that the lines are parallel, so the given system has no solution. It can be
proved algebraically simply by multiplying the first equation of the system
by 2 to get a system of the form
2x1 + 2x2 = 2
2x1 + 2x2 = 3,
Second, the two lines may not be parallel, and they may meet at exactly
one point, so in this case the system has exactly one solution. For example,
consider the system
x1 − x2 = −1
3x1 − x2 = 3.
From the graphs (Figure 1.1(b)) of these two equations we can see that
the lines intersect at exactly one point, namely, (2, 3), and so the system
has exactly one solution, x1 = 2, x2 = 3. To show this algebraically, if we
substitute x2 = x1 + 1 in the second equation, we have 3x1 − x1 − 1 = 3,
or x1 = 2, and using this value of x in x2 = x1 + 1 gives x2 = 3.
6 Applied Linear Algebra and Optimization using MATLAB
Finally, the two lines may actually be the same line, and so in this case,
every point on the lines gives a solution to the system and therefore there
are infinitely many solutions. For example, consider the system
x1 + x2 = 1
2x1 + 2x2 = 2.
Here, both equations have the same line for their graph (Figure 1.1(c)).
So this system has infinitely many solutions because any point on this line
gives a solution to this system, since any solution of the first equation is
also a solution of the second equation. For example, if we set x2 = x1 − 1,
and choose x1 = 0, x2 = 1, x1 = 1, x2 = 0, and so on. •
Note that a system of equations with no solution is said to be an incon-
sistent system and if it has at least one solution, it is said to be a consistent
system.
Matrices and Linear Systems 7
The system of linear equations (1.3) can be written as the single matrix
equation
a11 a12 · · · a1n x1 b1
a21 a22 · · · a2n x2 b2
.. .. = .. . (1.4)
.. .. ..
. . . . . .
an1 an2 · · · ann xn bn
If we compute the product of the two matrices on the left-hand side of
(1.9), we have
a11 x1 + a12 x2 + · · · + a1n xn b1
a21 x1 + a22 x2 + · · · + a2n xn b2
= . (1.5)
.. .. .. .. ..
. . . . .
an1 x1 + an2 x2 + · · · + ann xn bn
But two matrices are equal if and only if their corresponding elements
are equal. Hence, the single matrix equation (1.9) is equivalent to the
system of the linear equations (1.3). If we define
a11 a12 · · · a1n x1 b1
a21 a22 · · · a2n x2 b2
A = .. , x = , b = .. ,
.. .. .. ..
. . . . . .
an1 an2 · · · ann xn bn
the coefficient matrix, the column matrix of unknowns, and the column
matrix of constants, respectively, and then the system (1.3) can be written
very compactly as
Ax = b, (1.6)
8 Applied Linear Algebra and Optimization using MATLAB
which is called the matrix form of the system of linear equations (1.3). The
column matrices x and b are called vectors.
If the right-hand sides of the equal signs of (1.6) are not zero, then the
linear system (1.6) is called a nonhomogeneous system, and we will find
that all the equations must be independent to obtain a unique solution.
If the constants b of (1.6) are added to the coefficient matrix A as a
column of elements in the position shown below
..
a11 a12 · · · a1n . b1
..
a21 a22 · · · a2n . b2
[A|b] = , (1.7)
.. .. .. .. .. ..
. . . . . .
..
an1 an2 · · · ann . bn
then the matrix [A|b] is called the augmented matrix of the system (1.6).
In many instances, it may be convenient to operate on the augmented ma-
trix instead of manipulating the equations. It is customary to put a bar
between the last two columns of the augmented matrix to remind us where
the last column came from. However, the bar is not absolutely necessary.
The coefficient and augmented matrices of a linear system will play key
roles in our methods of solving linear systems.
>> A = [1 2 3; 4 5 6; 7 8 9];
>> b = [10; 11; 12];
>> Aug = [A b]
Aug =
1 2 3 10
4 5 6 11
7 8 9 12
Also,
Matrices and Linear Systems 9
The system of linear equations (1.8) can be written as the single matrix
equation
a11 a12 · · · a1n x1 0
a21 a22 · · · a2n x2 0
.. .. = .. . (1.9)
.. .. ..
. . . . . .
an1 an2 · · · ann xn 0
It can also be written in more compact form as
Ax = 0, (1.10)
where
a11 a12 · · · a1n x1 0
a21 a22 · · · a2n x2 0
A= .. , x= , 0= .
.. .. .. .. ..
. . . . . .
an1 an2 · · · ann xn 0
system, there are only two possibilities: either the zero solution is the only
solution, or there are infinitely many solutions (called nontrivial solutions).
Of course, it is usually nontrivial solutions that are of interest in physical
problems. A nontrivial solution to the homogeneous system can occur with
certain conditions on the coefficient matrix A, which we will discuss later.
The numbers a11 , a12 , . . . , amn that make up the array are called the ele-
ments of the matrix. The first subscript for the element denotes the row
and the second denotes the column in which the element appears. The
elements of a matrix may take many forms. They could be all numbers
(real or complex), or variables, or functions, or integrals, or derivatives, or
even matrices themselves.
The order or size of a matrix is specified by the number of rows (m)
and column (n); thus, the matrix A in (1.11) is of order m by n, usually
written as m × n.
A vector can be considered a special case of a matrix having only one
row or one column. A row vector containing n elements is a 1 × n matrix,
called a row matrix, and a column vector of n elements is an n × 1 matrix,
called a (column matrix). A matrix of order 1 × 1 is called a scalar.
Matrices and Linear Systems 11
Two matrices A = (aij ) and B = (bij ) are equal if they are the same size
and the corresponding elements in A and B are equal, i.e.,
Then
1 2 4 1 5 3
+ = = C.
3 4 5 2 8 6
•
>> A = [1 2; 3 4];
>> B = [4 1; 5 2];
>> C = A + B
C=
5 3
8 6
Then
1 2 4 1 4 + 10 1 + 4 14 5
= = = C.
3 4 5 2 12 + 20 3 + 8 32 11
then
2 3 1 7
AB = while BA = .
−2 2 −1 3
Thus, AB 6= BA. •
>> A = [1 2; −1 3];
>> B = [2 1; 0 1];
>> C = A ∗ B
C=
2 3
−2 2
>> u = [1 2 3 4];
>> v = [5 3 0 2];
>> x = u. ∗ v
x=
5608
>> y = u./v;
y=
0.2000 0.6667 Inf 2.0000
>> A = [1 2 3; 4 5 6; 7 8 9];
>> B = [9 8 7; 6 5 4; 3 2 1];
>> C = A. ∗ B
C=
9 16 21
24 25 24
21 16 9
>> A = [1 2 3; 4 5 6; 7 8 9];
>> D = A.ˆ 2
D=
1 4 9
16 25 36
49 64 81
Matrices and Linear Systems 15
A matrix A which has the same number of rows m and columns n, i.e.,
m = n, defined as
are square matrices because both have the same number of rows and columns.
•
16 Applied Linear Algebra and Optimization using MATLAB
In MATLAB, identity matrices are created with the eye function, which
can take either one or two input arguments:
>> I = eye(n)
>> I = eye(m, n)
>> A = [1 2 3; 4 5 6; 7 8 9]
>> B = A0
B=
1.0000 4.0000 7.0000
2.0000 5.0000 8.0000
3.0000 6.0000 9.0000
Note that
1. (AT )T = A,
AB = BA = In .
Then the matrix B is called the inverse of A and is denoted by A−1 . For
example, let
−1 32
2 3
A= and B= .
2 2 1 1
Then we have
AB = BA = I2 ,
>> A = [21 0 0; −1 2 − 1 0; 0 − 1 2 − 1; 0 0 − 1 2]
>> Ainv = IN V M AT (A)
Ainv =
0.8000 0.6000 0.4000 0.2000
0.6000 1.2000 0.8000 0.4000
0.4000 0.8000 1.2000 0.6000
0.2000 0.4000 0.6000 0.8000
Matrices and Linear Systems 19
Program 1.1
MATLAB m-file for Finding the Inverse of a Matrix
function [Ainv]=INVMAT(A)
[n,n]=size(A); I=zeros(n,n);
for i=1:n; I(i,i)=1; end
m(1:n,1:n)=A; m(1 : n, n + 1 : 2 ∗ n) = I;
for i=1:n; m(i, 1 : 2 ∗ n) = m(i, 1 : 2 ∗ n)/m(i, i);
for k=1:n; if i˜ =k
m(k, 1 : 2∗n) = m(k, 1 : 2∗n)−m(k, i)∗m(i, 1 : 2∗n);
end; end; end
invrs = m(1 : n, n + 1 : 2 ∗ n);
>> I = Ainv ∗ A;
>> f ormat short e
>> disp(I)
I=
1.0000e + 00 −1.1102e − 16 0 0
0 1.0000e + 00 0 0
0 0 1.0000e + 00 2.2204e − 16
0 0 0 1.0000e + 00
The values of I(2, 1), and I(3, 4) are very small, but nonzero, due to
round-off errors in the computation of Ainv and I. It is often preferable to
use rational numbers rather than decimal numbers. The function frac(x)
returns the rational approximation to x, or we can use the other MATLAB
command as follows:
3. Its product with another invertible matrix is invertible, and the in-
verse of the product is the product of the inverses in the reverse order.
If A and B are invertible matrices of the same size, then AB is in-
vertible and (AB)−1 = B −1 A−1 .
It is a square matrix having all elements equal to zero except those on the
main diagonal, i.e.,
aij = 0, if i 6= j
A = (aij ) =
aij 6= 0, if i = j.
Note that all diagonal matrices are invertible if all diagonal entries are
nonzero. •
Matrices and Linear Systems 21
>> B = [2 − 4 1; 6 10 − 3; 0 5 8]
>> M = diag(B)
M=
2
10
8
It is a square matrix which has zero elements below and to the left of the
main diagonal. The diagonal as well as the above diagonal elements can
take on any value, i.e.,
1 2 3
U = 0 4 5 .
0 0 6
>> A = [1 2 3; 4 5 6; 7 8 9];
>> U = triu(A)
U=
1 2 3
0 4 5
0 0 6
We can also create a strictly upper-triangular matrix, i.e., an upper-
triangular matrix with zero diagonals, from a given matrix A by using the
MATLAB built-in function triu(A,I) as follows:
>> A = [1 2 3; 4 5 6; 7 8 9];
>> U = triu(A, I)
U=
0 2 3
0 0 5
0 0 0
Matrices and Linear Systems 23
It is a square matrix which has zero elements above and to the right of the
main diagonal, and the rest of the elements can take on any value, i.e.,
Note that all the triangular matrices (upper or lower) with nonzero
diagonal entries are invertible.
>> A = [1 : 4; 5 : 8; 9 : 12]
%A is not symmetric
>> B = A0 ∗ A
B=
107 122 137 152
122 140 158 176
137 158 179 200
152 176 200 224
>> C = A ∗ A0
C=
30 70 110
70 174 278
110 278 446
Example 1.1 Find all the values of a, b, and c for which the following
matrix is symmetric:
4 a+b+c 0
A = −1 3 b − c .
−a + 2b − 2c 1 b − 2c
Solution. If the given matrix is symmetric, then A = AT , i.e.,
4 a+b+c 0
A = −1 3 b−c
−a + 2b − 2c 1 b − 2c
4 −1 −a + 2b − 2c
= a+b+c 3 1 = AT ,
0 b − c b − 2c
Matrices and Linear Systems 25
0 = −a + 2b − 2c
−1 = a + b + c
1 = b − c.
a = 2, b = −1, c = −2,
and using these values, we have the given matrix of the form
4 −1 0
A= −1 3 1 .
0 1 3
•
Theorem 1.3 If A and B are symmetric matrices of the same size, and
if k is any scalar, then:
1. AT is also symmetric;
If for a matrix A, the aij = −aji for a i 6= j and the main diagonal
elements are not all zero, then the matrix A is called a skew matrix. If
all the elements on the main diagonal of a skew matrix are zero, then the
matrix is called skew symmetric, i.e.,
Any square matrix may be split into the sum of a symmetric and a skew
symmetric matrix. Thus,
1 1
A = (A + AT ) + (A − AT ),
2 2
1
where 2 (A+A ) is a symmetric matrix and 12 (A−AT ) is a skew symmetric
T
The partitioning lines must always extend entirely through the matrix
as in the above example. If the submatrices of A are denoted by the symbols
A11 , A12 , A21 , and A22 so that
a11 a12 a13 a14 a15
A11 = a21 a22 a23 , A12 = a24 a25 ,
a31 a32 a33 a34 a35
A21 = a41 a42 a43 , A22 = a44 a45 ,
then the original matrix can be written in the form
A11 A12
A= .
A21 A22
A partitioned matrix may be transposed by appropriate transposition
and rearrangement of the submatrices. For example, it can be seen by
inspection that the transpose of the matrix A is
T
A11 AT12
AT = .
T T
A21 A22
Partitioned matrices such as the one given above can be added, sub-
tracted, and multiplied provided that the partitioning is performed in an
appropriate manner. For the addition and subtraction of two matrices, it is
necessary that both matrices be partitioned in exactly the same way. Thus,
a partitioned matrix B of order 4 × 5 (compare with matrix A above) will
be conformable for addition with A only if it is partitioned as follows:
.
b11 b12 b13 .. b14 b15
b21 b22 b23 ... b24 b25
B = b31 b32 b33 ... b34 b35 .
..
··· ··· ··· . ··· ···
.
b41 b42 b43 .. b44 b45
28 Applied Linear Algebra and Optimization using MATLAB
in which B11 , B12 , B21 , and B22 represent the corresponding submatrices.
In order to add A and B and obtain a sum C, it is necessary according to
the rules for addition of matrices that the following represent the sum:
A11 + B11 A12 + B12 C11 C12
A+B = = = C.
A21 + B21 A22 + B22 C21 C22
Note that like A and B, the sum matrix C will also have the same par-
titions.
Then, when forming the product AD according to the usual rules for
matrix multiplication, the following result is obtained:
A11 A12 D11 D12
M = AD =
A21 A22 D21 D22
A11 D11 + A12 D21 A11 D12 + A12 D22
=
A21 D11 + A22 D21 A21 D12 + A22 D22
M11 M12
= .
M21 M22
same way as the rows of the second partitioned matrix. It does not matter
how the rows of the first partitioned matrix and the columns of the second
partitioned matrix are partitioned. •
0 0 7 8
is a tridiagonal matrix.
A permutation matrix P has only 0s and 1s and there is exactly one in each
row and column of P . For example, the following matrices are permutation
matrices:
0 1 0 0
1 0 0 1 0 0 0
P = 0 0 1 , P = 0 0 1 0 .
0 1 0
0 0 0 1
The product P A has the same rows as A but in a different order (permuted),
while AP is just A with the columns permuted. •
2. The first entry from the left of a nonzero row is 1. This entry is
called the leading one of its row.
Matrices and Linear Systems 31
3. For each nonzero row, the leading one appears to the right and below
any leading ones in preceding rows.
Note that, in particular, in any column containing a leading one, all
entries below the leading one are zero. For example, the following matrices
are in row echelon form:
0 1 0 1
1 2 1 1 0 2 1 2 3 4
0 1 3 , 0 1 2 , 0 0 1 2 , 0 0 1 0 .
0 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0
0x1 + 0x2 = 0,
Similarly, the linear system that corresponds to the second above matrix
is
x1 = 2
x2 = 2
0 = 1.
The third equation of this system shows that
0x1 + 0x2 = 1,
32 Applied Linear Algebra and Optimization using MATLAB
which is not possible for any choices of x1 and x2 . Hence, the system has
no solution.
Finally, the linear system that corresponds to the third above matrix is
x1 + 2x2 + 3x3 = 4
x3 = 2
0x1 + 0x2 + 0x3 = 0,
and by backward substitution (without using the third equation of the sys-
tem), we get
2. The first entry from the left of a nonzero row is 1. This entry is
called the leading one of its row.
3. For each nonzero row, the leading one appears to the right and below
any leading ones in preceding rows.
For example, the following matrices are in reduced row echelon form:
1 0 1 1 0 0 2 1 4 5 0 1 1 0 0 0
0 1 2 , 0 1 0 4 , 0 0 0 1 , 0 0 1 0 2 ,
0 0 0 0 0 1 6 0 0 0 0 0 0 0 1 1
and the following matrices are not in reduced row echelon form:
1 3 0 2 1 3 0 2 1 0 0 3 1 0 2 0 0
0 0 0 0 , 0 0 5 4 , 0 0 1 2 , 0 1 1 0 2 .
0 0 1 4 0 0 0 1 0 1 0 6 0 2 0 1 1
There are usually many sequences of row operations that can be used to
transform a given matrix to reduced row echelon form—they all, however,
lead to the same reduced row echelon form. In the following, we shall
discuss how to transform a given matrix in reduced row echelon form.
Theorem 1.4 Every matrix can be brought to reduced row echelon form
by a series of elementary row operations. •
Using the finite sequence of elementary row operations, we get the matrix
of the form
1 −3 0 0 1
R1 = 0 0 1 −1 1 ,
0 0 0 1 0
Matrices and Linear Systems 35
which is in row echelon form. If we continue with the matrix R1 and make
all elements above the leading one equal to zero, we obtain
1 −3 0 0 1
R2 = 0 0 1 0 1 ,
0 0 0 1 0
>> A = [1 − 3 0 0 1; 2 − 6 − 1 1 1; 3 − 9 2 − 1 5];
>> B = rref (A)
B=
1 −3 0 0 1
0 0 1 0 1
0 0 0 1 0
so R1 is row equivalent to A.
Interchanging row 2 and row 3 of the matrix R1 gives the matrix of the
form
1 3 6 5
R2 = 2 −7 −3 −1 ,
2 1 4 3
so R2 is row equivalent to R1 .
Theorem 1.5
1. Every matrix is row equivalent to itself.
2. If a matrix A is row equivalent to a matrix B, then B is row equivalent
to A.
3. If a matrix A is row equivalent to a matrix B and B is row equivalent
to a matrix C, then A is row equivalent to C. •
Example 1.5 Use elementary row operations on matrices to solve the lin-
ear system
− x2 + x3 = 1
x1 − x2 − x3 = 1
−x1 + 3x3 = −2.
Matrices and Linear Systems 37
x1 = −4
x2 = −3
x3 = −2,
and we get the solution [−4, −3, −2]T of the given linear system. •
1. det(A) = a11 , if n = 1.
For example, if
4 2 6 3
A= and B = ,
−3 7 2 5
then
>> A = [2 2; 6 7];
>> B = det(A)
B=
2.0000
Another way to find the determinants of only 2 × 2 and 3 × 3 matrices
can be found easily and quickly using diagonals (or direct evaluation). For
a 2 × 2 matrix, the determinant can be obtained by forming the product of
the entries on the line from left to right and subtracting from this number
the product of the entries on the line from right to left. For a matrix of
size 3 × 3, the diagonals of an array consisting of the matrix with the first
two columns added to the right are used. Then the determinant can be
obtained by forming the sum of the products of the entries on the lines
from left to right, and subtracting from this number the products of the
entries on the lines from right to left, as shown in Figure (1.2).
then the minor M11 will be obtained by deleting the first row and the first
column of the given matrix A, i.e.,
3 2
M11 = = (3)(4) − (−2)(2) = 12 + 4 = 16.
−2 4
Similarly, we can find the other possible minors of the given matrix as
Matrices and Linear Systems 41
follows:
5 2
M12 =
4 4
=
20 − 16 = 4
5 3
M13 = = −10 − 12 = −22
4 −2
3 −1
M21 = = 12 − 2 = 10
−2 4
2 −1
M22 =
4
= 8+4 = 12
4
2 3
M23 =
4 −2
=
−4 − 12 = −16
3 −1
M31 =
3
= 6+3 = 9
2
2 −1
M32 =
5
= 4+5 = 9
2
2 3
M33 = 5 3 =
6 − 15 = −9,
>> A = [2 3 − 1; 5 3 2; 4 − 2 4];
>> Cof A = cof actor(A, 1, 1);
>> Cof A = cof actor(A, 1, 2);
>> Cof A = cof actor(A, 1, 3);
>> Cof A = cof actor(A, 2, 1);
>> Cof A = cof actor(A, 2, 2);
>> Cof A = cof actor(A, 2, 3);
>> Cof A = cof actor(A, 3, 1);
>> Cof A = cof actor(A, 3, 2);
>> Cof A = cof actor(A, 3, 3);
Matrices and Linear Systems 43
Program 1.2
MATLAB m-file for Finding Minors and Cofactors
of a Matrix
function CofA = cofactor(A,i,j)
[m,n] = size(A);
if m ˜ = n error(Matrix must be square) end
A1 = A([1:i-1,i+1:n],[1:j-1,j+1:n]);
Minor = det(A1);
CofA = (-1)ˆ (i+j)*det(Minor);
where the summation is on i for any fixed value of the jth column (1 ≤ j ≤
n), or on j for any fixed value of the ith row (1 ≤ i ≤ n), and Aij is the
cofactor of element aij . •
Example 1.6 Find the minors and cofactors of the matrix A and use them
to evaluate the determinant of the matrix
3 1 −4
A= 2 5 6 .
1 4 8
44 Applied Linear Algebra and Optimization using MATLAB
From these values of the minors, we can calculate the cofactors of the
elements of the given matrix as follows:
Now by using the cofactor expansion along the first row, we can find
the determinant of the matrix as follows:
>> A = [3 1 − 4; 2 5 6; 1 4 8];
>> DetA = Cof F exp(A);
Matrices and Linear Systems 45
Program 1.3
MATLAB m-file for Finding the Determinant of a
Matrix by Cofactor Expansion
function DetA = CofFexp(A)
[m,n] = size(A);
if m ˜ = n error (Matrix must be square) end
a = A(1,:);c = [ ];
for i=1:n
c1i = cofactor(A,1,i);
c = [c;c1i]; end
DetA = a*c;;
which is called the cofactor expansion along the ith row, and also as
and is called the cofactor expansion along the jth column. This is called
the Laplace expansion theorem. •
Note that the cofactor and minor of an element aij differs only in sign,
i.e., Aij = ±Mij . A quick way for determining whether to use the + or −
is to use the fact that the sign relating Aij and Mij is in the ith row and
jth column of the checkerboard array
46 Applied Linear Algebra and Optimization using MATLAB
+ − + − + ···
− + − + − ···
+ − + − + ··· .
− + − + − ···
.. .. .. .. .. ..
. . . . . .
For example, A11 = M11 , A21 = −M21 , A12 = −M12 , A22 = M22 , and so on.
If A is any n × n matrix and Aij is the cofactor of aij , then the matrix
A11 A12 · · · A1n
A21 A22 · · · A2n
.. .. ..
. . ··· .
An1 An2 · · · Ann
is called the matrix of the cofactor from A. For example, the cofactor of
the matrix
3 2 −1
A= 1 6 3
2 −4 0
can be calculated as follows:
Thus,
The following are special properties, which will be helpful in reducing the
amount of work involved in evaluating determinants.
Let A be an n × n matrix:
det(B) = 8 = det(A).
then
det(A) = (3)(4)(5) = 60.
11. The determinant of the kth power of a matrix A is equal to the kth
power of the determinant of the matrix A, i.e., det(Ak ) = (det(A))k .
For example, if
2 −2 0
A= 2 3 −1 ,
1 0 1
then det(A) = 12, and for the matrix
−18 −30 12
B = A3 = 24 −3 −9 ,
3 −12 3
obtained by taking the cubic power of the matrix A, we have
Example 1.8 Find all the values of α for which det(A) = 0, where
α−3 1 0
A= 0 α − 1 1 .
0 2 α
Solution. We find the determinant of the given matrix by using the co-
factor expansion along the first row, so we compute
|A| = 0
(α − 3)(α + 1)(α − 2) = 0,
which gives
α = −1, α = 2, α = 3,
the required values of α for which det(A) = 0. •
Solution. Since
4α α
= 4α(α + 1) − α,
1 α+1
which is equivalent to
4α α
= 4α2 + 3α.
1 α+1
Matrices and Linear Systems 53
Also,
3 −1 0
0 α
−2 = 3[α(α+1)+6]−(−1)[(0)(α+1)+2]+0[(0)(3)−(−1)(α)],
−1 3 α+1
Given that
3 −1 0
4α α
1 α+1
= 0 α
−2 ,
−1 3 α+1
we get
4α2 + 3α = 3α2 + 3α + 16.
Simplifying this quadratic polynomial, we have
α2 = 16 or α2 − 16 = 0,
which gives
α = −4 and α = 4,
the required values of α. •
if
a b c
det d e f = 4.
g h i
54 Applied Linear Algebra and Optimization using MATLAB
or
a b c
|A| = (−5)(−1)(2)(−1) d e f .
g h i
Since it is given that
a b c
d e f = 4,
g h i
we have
|A| = (−5)(−1)(2)(−1)(4) = −40,
the required determinant of the given matrix. •
One can easily transform the given determinant into upper-triangular form
by using the following row operations:
1. Add a multiple of one row to another row, and this will not affect
the determinant.
Matrices and Linear Systems 55
Now to create the zeros below the main diagonal, column by column, we
do as follows:
Replace the second row of the determinant with the sum of itself and (−6)
times the first row of the determinant and then replace the third row of
the determinant with the sum of itself and (3) times the first row of the
determinant, which gives
1 2 3
(3) 0 −10 −25 .
0 7 8
1
Multiplying row 2 of the determinant by − 10 gives
1 2 3
(3)(−10) 0 1 −5/2 .
0 7 8
Replacing the third row of the determinant with the sum of itself and (−7)
times the second row of the determinant, we obtain
56 Applied Linear Algebra and Optimization using MATLAB
1 2 3
(3)(−10) 0 1 −5/2
= (3)(−10)(1)(1)(−19/2) = 285,
0 0 −19/2
Example 1.12 For what values of α does the following matrix have an
inverse?
1 0 α
A= 2 2 1
0 2α 1
Solution. We find the determinant of the given matrix by using cofactor
expansion along the first row as follows:
which is equal to
2 1
C11 = (−1)1+1 M11 = M11 = = 2 − 2α
2α 1
2 2
C13 = (−1)1+3 M13 = M13 = − = 4α.
0 2α
Thus,
From Theorem 1.9 we know that the matrix has an inverse if det(A) 6= 0,
so
|A| = 2 − 2α + 4α2 6= 0,
Example 1.13 Use the adjoint method to compute the inverse of the fol-
lowing matrix:
1 2 −1
A = 2 −1 1 .
1 2 2
Also, find the inverse and determinant of the adjoint matrix.
−1 1 2 −1
= −4, C12 = − 2 1
C11 = + = −3, C13 = + = 5,
2 2 1 2 1 2
2 −1 −1
= −6, C22 = + 1
1 2
C21 = − = 3, C23 = − = 0,
2 2 1 2 1 2
2 −1 1 −1 1 2
C31 = + = 1, C32 = − = −3, C33 = + = −5.
−1 1 2 1 2 −1
Thus, the cofactor matrix has the form
−4 −3 5
−6 3 0 ,
1 −3 −5
and the adjoint is the transpose of the cofactor matrix
T
−4 −3 5 −4 −6 1
adj(A) = −6 3 0 = −3 3 −3 .
1 −3 −5 5 0 −5
To get the adjoint of the matrix of Example 1.13, we use the MATLAB
Command Window as follows:
>> A = [1 2 − 1; 2 − 1 1; 1 2 2];
>> AdjA = Adjoint(A);
Program 1.4
MATLAB m-file for Finding the Adjoint of a
Matrix Function AdjA = Adjoint(A)
[m,n] = size(A);
if m ˜ = n error(‘Matrix must be square’) end
A1 = [ ];
for i = 1:n
for j=1:n
A1 = [A1;cofactor(A,i,j)];end;end
AdjA = reshape(A1,n,n);
Matrices and Linear Systems 59
Then by using Theorem 1.9 we can have the inverse of the matrix as
follows:
−4 −6 1 4/15 2/5 −1/15
Adj(A) 1
A−1 = = − −3 3 −3 = 1/5 −1/5 1/5 .
det(A) 15
5 0 −5 −1/3 0 1/3
Using Theorem 1.9 we can compute the inverse of the adjoint matrix as:
−1/15 −2/15 1/15
A
(adj(A))−1 = = −2/15 1/15 −1/15 ,
det(A)
−1/15 −2/15 −2/15
by using the adjoint and the determinant of the matrix in the MATLAB
Command Window as:
>> A = [1 − 1 1 2; 1 0 1 3; 0 0 2 4; 1 1 − 1 1];
The cofactors Aij of elements of the given matrix A can also be found
directly by using the MATLAB Command Window as follows:
60 Applied Linear Algebra and Optimization using MATLAB
>> B =
[A11 A12 A13 A14;
A21 A22 A23 A24;
A31 A32 A33 A34;
A41 A42 A43 A44]
which gives
B=
−2 −4 −4 2
6 6 8 −4
−3 −2 −3 2
−2 −2 −4 2
>> adjA = B 0
−2 6 −3 −2
−4 6 −2 −2
−4 8 −3 −4
2 −4 2 2
The determinant of the matrix can be obtained as:
>> det(A)
ans =
2
The inverse of A is the adjoint matrix divided by the determinant of A.
>> inv(A)
det(A2 B −1 AT B 3 ) = 432.
The system of linear equations (1.14) can be written as the single matrix
equation
a11 a12 · · · a1n x1 0
a21 a22 · · · a2n x2 0
.. .. = .. . (1.15)
.. .. ..
. . . . . .
am1 am2 · · · amn xn 0
If we compute the product of the two matrices on the left-hand side of
(1.15), we have
a11 x1 + a12 x2 + · · · + a1n xn 0
a21 x1 + a22 x2 + · · · + a2n xn 0
= . (1.16)
.. .. .. .. ..
. . . . .
am1 x1 + am2 x2 + · · · + amn xn 0
But the two matrices are equal if and only if their corresponding elements
are equal. Hence, the single matrix equation (1.15) is equivalent to the
system of the linear equations (1.14). If we define
a11 a12 · · · a1n x1 0
a21 a22 · · · a2n x2 0
A = .. .. , x = .. , b = .. ,
.. ..
. . . . . .
am1 am2 · · · amn xn 0
the coefficient matrix, the column matrix of unknowns, and the column
matrix of constants, respectively, then the system (1.14) can be written
very compactly as
Ax = b, (1.17)
which is called the matrix form of the homogeneous system. •
[A|0].
64 Applied Linear Algebra and Optimization using MATLAB
x1 + x2 + 2x3 = 0
2x1 + 3x2 + 4x3 = 0
3x1 + 4x2 + 7x3 = 0.
Next, using the elementary row operations: row3 – row2 and row1 – row2,
we get
..
1 0 2 . 0
≡ 0 1 0 ... 0
.
..
0 0 1 . 0
Thus,
x1 = 0, x2 = 0, x3 = 0
x1 + 2x2 + x3 = 0
2x1 − 3x2 + 4x3 = 0.
. !
1 2 1 .. 0
∼ . .
0 1 − 72 .. 0
11
x1 + 0x2 + x3 = 0
7
2
0x1 + x2 − x3 = 0
7
and from it, we get
2
x2 = x3 .
7
Taking x3 = t, for t ∈ R and t 6= 0, we get the nontrivial solution
11 2
[x1 , x2 , x3 ]T = [ t, t, t]T .
7 7
Thus, the given system has infinitely many solutions, and this is to be
expected because the given system has three unknowns and only two equa-
tions. •
Example 1.17 For what values of α does the homogeneous linear system
(α − 2)x1 + x2 = 0
x1 + (α − 2)x2 = 0
.. !
1 (α − 2) . 0
∼ .. .
0 1 − (α − 2)2 . 0
Using backward substitution, we obtain
x1 + (α − 2)x2 = 0
0x1 + 1 − (α − 2)2 x2 = 0.
Notice that if x2 = 0, then x1 = 0, and the given system has a trivial
6 0. This implies that
solution, so let x2 =
1 − (α − 2)2 = 0
1 − α2 + 4α − 4 = 0
α2 − 4α + 3 = 0
(α − 3)(α − 1) = 0,
which gives
α = 1 and α = 3.
Notice that for these values of α, the given set of equations are identical,
i.e.,
(for α = 1)
−x1 + x2 = 0
x1 − x2 = 0,
and (for α = 3)
x1 + x2 = 0
x1 + x2 = 0.
Thus, the given system has nontrivial solutions (infinitely many solu-
tions) for α = 1 and α = 3. •
The following basic theorems on the solvability of linear systems are
proved in linear algebra.
68 Applied Linear Algebra and Optimization using MATLAB
A−1 Ax = A−1 b
Ix = A−1 b
or
x = A−1 b. (1.18)
If A is a square invertible matrix, there exists a sequence of elementary
row operations that carry A to the identity matrix I of the same size, i.e.,
A −→ I. This same sequence of row operations carries I to A−1 , i.e.,
I −→ A−1 . This can also be written as
[A|I] −→ [I|A−1 ].
Example 1.18 Use the matrix inversion method to find the solution of
the following linear system:
x1 + 2x2 = 1
−2x1 + x2 + 2x3 = 1
−x1 + x2 + x3 = 1.
Matrices and Linear Systems 69
Multiply the first row by −2 and −1 and then, subtracting the results
from the second and third rows, respectively, we get
.
1 2 0 .. 1 0 0
∼
..
.
0 5 2 . 2 1 0
..
0 3 1 . 1 0 1
Multiplying the second row by 2 and 3 and then subtracting the results
from the first and third rows, respectively, we get
.
1 0 −4/5 .. 1/5 −2/5 0
.
∼ 0 1 2/5 .. 2/5 1/5 .
0
..
0 0 −1/5 . −1/5 −3/5 1
.
1 0 −4/5 .. 1/5 −2/5 0
∼ ..
.
0 1 2/5 . 2/5 1/5 0
..
0 0 1 . 1 3 −5
Multiplying the third row by 25 and − 45 and then subtracting the results
from the second and first rows, respectively, we get
..
1 0 0 . 1 2 −4
.
∼ 0 1 0 .. 0 −1
.
2
.
0 0 1 .. 1 3 −5
Thus, the inverse of the given matrix is
1 2 −4
A−1 = 0 −1 2 ,
1 3 −5
and the unique solution of the system can be computed as
1 2 −4 1 −1
x = A−1 b = 0 −1 2 1 = 1 ,
1 3 −5 1 −1
i.e.,
x1 = −1, x2 = 1, x3 = −1,
the solution of the given system by the matrix inversion method. •
Thus, when the matrix inverse A−1 of the coefficient matrix A is com-
puted, the solution vector x of the system (1.6) is simply the product of
inverse matrix A−1 and the right-hand side vector b.
>> A = [1 2 0; −2 1 2; −1 1 1];
>> b = [1; 1; 1];
>> x = A \ b
x=
1.0000
1.0000
−1.0000
Not all matrices have inverses. Singular matrices don’t have inverses
and thus the corresponding systems of equations do not have unique solu-
tions. The inverse of a matrix can also be computed by using the following
numerical methods for linear systems: Gauss-elimination method, Gauss–
Jordan method, and LU decomposition method. But the best and simplest
method for finding the inverse of a matrix is to perform the Gauss–Jordan
method on the augmented matrix with an identity matrix of the same size.
E3 E2 E1 A = I,
and so
A = (E3 E2 E1 )−1 .
This means that
0 1 1 0 1 1
A= E1−1 E2−1 E3−1 = .
1 0 2 1 0 1
•
74 Applied Linear Algebra and Optimization using MATLAB
4. It has n pivots. •
In the following, we will discuss the direct methods for solving the linear
systems.
In a similar way, one can use Cramer’s rule for a set of n linear equations
as follows:
76 Applied Linear Algebra and Optimization using MATLAB
|Ai |
xi = , i = 1, 2, 3, . . . , n, (1.19)
|A|
i.e., the solution for any one of the unknown xi in a set of simultaneous
equations is equal to the ratio of two determinants; the determinant in
the denominator is the determinant of the coefficient matrix A, while the
determinant in the numerator is the same determinant with the ith column
replaced by the elements from the right-hand sides of the equation.
5x1 + x3 + 2x4 = 3
x1 + x2 + 3x3 + x4 = 5
x1 + x2 + 2x4 = 1
x1 + x2 + x3 + x4 = −1.
gives
5 0 1 2 3
1 1 3 1 5
A=
1
and b=
1 .
1 0 2
1 1 1 1 −1
The determinant of the matrix A can be calculated by using cofactor
expansion as follows:
5 0 1 2
1 1 3 1
|A| =
1 1 0 2
1 1 1 1
= a11 c11 + a12 c12 + a13 c13 + a14 c14 = 5(2) + 0(−2) + 1(0) + 2(0) = 10 6= 0,
Matrices and Linear Systems 77
which shows that the given matrix A is nonsingular. Then the matrices
A1 , A2 , A3 , and A4 can be computed as
3 0 1 2 5 3 1 2
5 1 3 1 1 5 3 1
A1 = 1 1 0 2
, A2 =
1
,
1 0 2
−1 1 1 1 1 −1 1 1
5 0 3 2 5 0 1 3
1 1 5 1 1 1 3 5
A3 = 1 1
, A4 = .
1 2 1 1 0 1
1 1 −1 1 1 1 1 −1
The determinant of the matrices A1 , A2 , A3 , and A4 can be computed
as follows:
|A2 | 70
x2 = = − = −7
|A| 10
|A3 | 30
x3 = = = 3
|A| 10
|A3 | 50
x4 = = = 5,
|A| 10
which is the required solution of the given system. •
3
system of n linear equations by Cramer’s rule will require N = (n + 1) n3
multiplications. Therefore, this rule is much less efficient for large values of
n and is at most never used for computational purposes. When the number
of equations is large (n > 4), other methods of solutions are more desirable.
Use MATLAB commands to find the solution of the above linear sys-
tem by Cramer’s rule as follows:
>> A = [5 0 1 2; 1 1 3 1; 1 1 0 2; 1 1 1 1];
>> b = [3; 5; 1; −1];
>> A1 = [b A(:, [2 : 4])];
>> x1 = det(A1)/det(A);
>> A2 = [A(:, 1) b A(:, [3 : 4])];
>> x2 = det(A2)/det(A);
>> A3 = [A(:, [1 : 2]) b A(:, 4)];
>> x3 = det(A3)/det(A);
>> A4 = [A(:, [1 : 3]) b];
>> x4 = det(A4)/det(A);
The m-file CRule.m and the following MATLAB commands can be used
to generate the solution of Example 1.21 as follows:
Matrices and Linear Systems 79
>> A = [5 0 1 2; 1 1 3 1; 1 1 0 2; 1 1 1 1];
>> b = [3; 5; 1; −1];
>> sol = CRule(A, b);
Program 1.5
MATLAB m-file for Cramer’s Rule for a Linear System
function sol=CRule(A,b)
[m, n] = size(A);
if m ˜ = n error(‘Matrix is not square.’); end
if det(A) == 0 error(‘Matrix is singular.’);end
for i = 1:n
B = A; B(:, i) = b;
sol(i) = det(B) / det(A);end
sol = sol’;
Forward Elimination
as the first pivotal equation with the first pivot element a11 . Then the first
equation times multiples mi1 = (ai1 /a11 ), i = 2, 3, . . . , n is subtracted from
the ith equation to eliminate the first variable x1 , producing an equivalent
system
an equivalent system
a11 x1 + a12 x2 + a13 x3 + · · · + a1n xn = b1
(1) (1) (1) (1)
a22 x2 + a23 x3 + · · · + a2n xn = b2
(2) (2) (2)
a33 x3 + · · · + a3n xn = b3 (1.24)
.. .. .. ..
. . . .
(2) (2) (2)
an3 x3 + · · · + ann xn = bn .
Backward Substitution
After the triangular set of equations has been obtained, the last equation of
system (1.26) yields the value of xn directly. The value is then substituted
into the equation next to the last one of the system (1.26) to obtain a value
of xn−1 , which is, in turn, used along with the value of xn in the second
82 Applied Linear Algebra and Optimization using MATLAB
The Gaussian elimination can be carried out by writing only the co-
efficients and the right-hand side terms in a matrix form, the augmented
matrix form. Indeed, this is exactly what a computer program for Gaus-
sian elimination does. Even for hand calculations, the augmented matrix
form is more convenient than writing all sets of equations. The augmented
matrix is formed as follows:
a11 a12 a13 · · · a1n | b1
a21 a22 a23 · · · a2n | b2
a31 a32 a33 · · · a3n | b3
. (1.28)
.. .. .. .. ..
. . . . . |
an1 an2 an3 · · · ann | bn
.
1 2 1 .. 2
0 1 1 ... −3 .
..
0 1 3 . 3
(1) (1)
Since a22 = 1 6= 0, we eliminate the entry in the a32 position by subtracting
1
the multiple m32 = = 1 of the second row from the third row to get
1
.
1 2 1 .. 2
0 1 1 ... −3 .
..
0 0 2 . 6
Obviously, the original set of equations has been transformed to an
upper-triangular form. Since all the diagonal elements of the obtaining
upper-triangular matrix are nonzero, the coefficient matrix of the given
system is nonsingular, and the given system has a unique solution. Now
expressing the set in algebraic form yields
x1 + 2x2 + x3 = 2
x2 + x3 = −3
2x3 = 6.
Now using backward substitution, we get
2x3 = 6 gives x3 = 3
x2 = −x3 − 3 = −(3) − 3 = −6 gives x2 = −6
x1 = 2 − 2x2 − x3 = 2 − 2(−6) − 3 gives x1 = 11,
which is the required solution of the given system. •
The above results can be obtained using MATLAB commands as fol-
lows:
>> B = [1 2 1 2; 2 5 3 1; 1 3 4 5];
%B = [A|b] = Augmented matrix
>> x = W P (B);
>> disp(x)
Matrices and Linear Systems 85
Program 1.6
MATLAB m-file for Simple Gaussian Elimination Method
function x=WP(B)
[n,t]=size(B); U=B;
for k=1:n-1; for i=k:n-1; m=U(i+1,k)/U(k,k);
for j=1:t; U(i+1,j)=U(i+1,j)-m*U(k,j);end;end end
i=n; x(i,1)=U(i,t)/U(i,i);
for i=n-1:-1:1; s=0;
for k=n:-1:i+1; s = s + U (i, k) ∗ x(k, 1); end
x(i,1)=(U(i,t)-s)/U(i,i); end; B; U; x; end
Example 1.23 Solve the following linear system using the simple Gaus-
sian elimination method:
x2 + x3 = 1
x1 + 2x2 + 2x3 = 1
2x1 + x2 + 2x3 = 3.
To solve this system, the simple Gaussian elimination method will fail
immediately because the element in the first row on the leading diagonal,
the pivot, is zero. Thus, it is impossible to divide that row by the pivot
86 Applied Linear Algebra and Optimization using MATLAB
x2 = 1 − x3 = 1 − 4 = −3
Example 1.24 Solve the following linear system using the simple Gaus-
sian elimination method:
x1 + x2 + x3 = 3
2x1 + 2x2 + 3x3 = 7
x1 + 2x2 + 3x3 = 6.
The first elimination step is to eliminate the elements a21 = 2 and a31 =
1 from the second and third rows by subtracting the multiples m21 = 21 = 2
and m31 = 11 = 1 of row 1 from row 2 and row 3, respectively, which gives
.
1 1 1 .. 3
0 0 1 ... 1 .
..
0 1 2 . 3
We finished the second column. So the third row of the equivalent upper-
triangular system is
0x1 + 0x2 + 0x3 = b − a. (1.31)
First, if (1.31) has no constraint on unknowns x1 , x2 , and x3 , then the
upper-triangular system represents only two nontrivial equations, namely,
2x1 − x2 + 3x3 = 1
4x2 − 4x3 = 2a − 2
in the three unknowns. As a result, one of the unknowns can be chosen
arbitrarily, say x3 = x∗3 , then x∗2 and x∗1 can be obtained by using backward
substitution:
1
x∗2 = 1/2a − 1/2 − x∗3 ; x∗1 = (1 + 1/2a − 1/2 − 4x∗3 ).
2
Hence,
1
x∗ = [ (1 + 1/2a − 1/2 − 4x∗3 ), 1/2a − 1/2 − x∗3 , x∗3 ]T
2
is an approximation solution of the given system for any value of x∗3 for
any real value of a. Hence, the given linear system is consistent (infinitely
many solutions).
Then using the multiple m32 = 42 = 2 of the second row from the third row,
we get
..
1 1 −2 . 0
..
.
0 2 1 . 0
..
0 0 −1 . 0
Obviously, the original set of equations has been transformed to an
upper-triangular form. Thus, the system has the unique solution [0, 0, 0]T ,
i.e., the system has only the trivial solution. •
Example 1.27 Find the value of k for which the following homogeneous
linear system has nontrivial solutions by using the simple Gaussian elimi-
nation method:
2x1 − 3x2 + 5x3 = 0
−2x1 + 6x2 − x3 = 0
4x1 − 9x2 + kx3 = 0.
Solution. The process begins with the augmented matrix form
.
2 −3 5 .. 0
−2 ..
,
6 −1 . 0
..
4 −9 k . 0
which gives
.
2 −3 5 .. 0
0 ..
.
3 4 . 0
..
0 −3 k − 10 . 0
−3
Also, by using the multiple m32 = 3
= −1, we get
..
2 −3 5 . 0
.
0
3 4 .. 0
.
.
0 0 k − 6 .. 0
k − 6 = 0, which gives k = 6.
Example 1.28 Use the simple Gaussian elimination method to find all
the values of α which make the following matrix singular:
1 −1 α
A= 2 2 1 .
0 α −1.5
92 Applied Linear Algebra and Optimization using MATLAB
Solution. Apply the forward elimination step of the simple Gaussian elim-
ination on the given matrix A and eliminate the element a21 by subtracting
from the second row the appropriate multiple of the first row. In this case,
the multiple is given as
1 −1 α
0 4 1 − 2α .
0 α −1.5
We finished the first elimination step. The second elimination step is
(1)
to eliminate element a32 = α by subtracting a multiple m32 = α4 of row 2
from row 3, which gives
1 −1 α
0 4 1 − 2α .
α(1 − 2α)
0 0 −1.5 −
4
To show that the given matrix is singular, we have to set the third diagonal
element equal to zero (by Theorem 1.19), i.e.,
α(1 − 2α)
−1.5 − = 0.
4
After simplifying, we obtain
2α2 − α − 6 = 0.
Solving the above quadratic equation, we get
3
α=− and α = 2,
2
which are the possible values of α, which make the given matrix singular.•
Example 1.29 Use the smallest positive integer value of α to find the
unique solution of the linear system Ax = [1, 6, −4]T by the simple Gaus-
sian elimination method, where
1 −1 α
A= 2 2 1 .
0 α −1.5
Matrices and Linear Systems 93
Solution. Since we know from Example 1.28 that the given matrix A is
singular when α = − 32 and α = 2, to find the unique solution we take the
smallest positive integer value α = 1 and consider the augmented matrix
as follows:
.
1 −1 1 .. 1
2 ..
.
2 1 . 6
..
0 1 −1.5 . −4
Applying the forward elimination step of the simple Gaussian elimina-
tion on the given matrix A and eliminating the element a21 by subtracting
from the second row the appropriate multiple m21 = 2 of the first row gives
..
1 −1 1 . 1
.
0
4 −1 .. .
4
.
0 1 −1.5 .. −4
(1)
The second elimination step is to eliminate element a32 = 1 by subtracting
a multiple m32 = 14 of row 2 from row 3, which gives
..
1 −1 1 . 1
..
.
0
4 −1 . 4
..
0 0 −5/4 . −5
x1 − x2 + x3 = 1
4x2 − x3 = 4
−5/4x3 = −5.
x3 = 4, x2 = 2, x1 = −1,
Note that the inverse of the nonsingular matrix A can be easily determined
by using the simple Gaussian elimination method. Here, we have to con-
sider the augmented matrix as a combination of the given matrix A and the
identity matrix I (the same size as A). To find the inverse matrix BA−1 ,
we must solve the linear system in which the jth column of the matrix B
is the solution of the linear system with the right-hand side the jth column
of the matrix I.
Example 1.30 Use the simple Gaussian elimination method to find the
inverse of the following matrix:
2 −1 3
A = 4 −1 6 .
2 −3 4
Solution. Suppose that the inverse A−1 = B of the given matrix exists
and let
2 −1 3 b11 b12 b13 1 0 0
AB = 4 −1 6 b21 b22 b23 = 0 1 0 = I.
2 −3 4 b31 b32 b33 0 0 1
which gives
b11 = 7, b21 = −2, b31 = −5.
Similarly, the solution of the second linear system
2 −1 3 b12 0
0 1 0 b22 = 1
0 0 1 b32 2
which gives
b12 = −5/2, b22 = 1, b32 = 2.
96 Applied Linear Algebra and Optimization using MATLAB
and it gives
b13 = −3/2, b23 = 0, b33 = 1.
Hence, the elements of the inverse matrix B are
5 3
7 − −
2 2
B = A−1 = −2 ,
1 0
−5 2 1
2. Check the first pivot element a11 6= 0, then move to the next step;
otherwise, interchange rows so that a11 6= 0.
3. Multiply row one by multiplier mi1 = ai1 /a11 and subtract to the ith
row for i = 2, 3, . . . , n.
4. Repeat steps 2 and 3 for the remaining pivots elements unless coeffi-
cient matrix A becomes upper-triangular matrix U .
bn−1
5. Use backward substitution to solve xn from the nth equation xn = n
ann
and solve the other (n − 1) unknown variables by using (1.27).
Matrices and Linear Systems 97
We finished the first elimination step. The second pivot is in the (2, 2)
position, but after eliminating the element below it, we find the triangular
form to be
1 2 4
0 −1 1 .
0 0 3
Since the number of pivots are three, the rank of the given matrix is 3. Note
that the original matrix is nonsingular since the rank of the 3 × 3 matrix
is 3. •
>> A = [1 2 4; 1 1 5; 1 1 6];
>> rank(A)
ans =
3
Note that:
rank(AB) ≤ min(rank(A), rank(B))
rank(A + B) ≤ rank(A) + rank(B)
rank(AAT ) = rank(A) = rank(AT A)
Although the rank of a matrix is very useful to categorize the behavior
of matrices and systems of equations, the rank of a matrix is usually not
computed. •
0.000100x1 + x2 = 1
x1 + x2 = 2,
which has the exact solution x = [1.00010, 0.99990]T . Now we solve this
system by simple Gaussian elimination. The first elimination step is to
eliminate the first variable x1 from the second equation by subtracting
Matrices and Linear Systems 99
multiple m21 = 10000 of the first equation from the second equation, which
gives
0.000100x1 + x2 = 1
− 10000x2 = −10000.
Partial Pivoting
Here, we develop an implementation of Gaussian elimination that utilizes
the pivoting strategy discussed above. In using Gaussian elimination by
partial pivoting (or row pivoting), the basic approach is to use the largest
(in absolute value) element on or below the diagonal in the column of
current interest as the pivotal element for elimination in the rest of that
column.
One immediate effect of this will be to force all the multiples used to be
not greater than 1 in absolute value. This will inhibit the growth of error in
the rest of the elimination phase and in subsequent backward substitution.
At stage k of forward elimination, it is necessary, therefore, to be able to
identify the largest element from |akk |, |ak+1,k |, . . . , |ank |, where these aik s
are the elements in the current partially triangularized coefficient matrix. If
this maximum occurs in row p, then the pth and kth rows of the augmented
matrix are interchanged and the elimination proceeds as usual. In solving
n linear equations, a total of N = n(n+1)/2 coefficients must be examined.
Example 1.32 Solve the following linear system using Gaussian elimina-
tion with partial pivoting:
x1 + x2 + x3 = 1
2x1 + 3x2 + 4x3 = 3
4x1 + 9x2 + 16x3 = 11.
Solution. For the first elimination step, since 4 is the largest absolute
coefficient of the first variable x1 , the first row and the third row are inter-
changed, which gives us
4x1 + 9x2 + 16x3 = 11
2x1 + 3x2 + 4x3 = 3
x1 + x2 + x3 = 1.
Eliminate the first variable x1 from the second and third rows by subtracting
the multiples m21 = 24 and m31 = 41 of row 1 from row 2 and row 3,
respectively, which gives
4x1 + 9x2 + 16x3 = 11
− 3/2x2 − 4x3 = −5/2
− 5/4x2 − x3 = −7/5.
Matrices and Linear Systems 101
For the second elimination step, − 32 is the largest absolute coefficient of the
second variable x2 , so eliminate the second variable x2 from the third row
by subtracting the multiple m32 = 56 of row 2 from row 3, which gives
x1 = 1, x2 = −1, x3 = 1,
The following MATLAB commands will give the same results we ob-
tained in Example 1.32 of the Gaussian elimination method with partial
pivoting:
>> B = [1 1 1 1; 2 3 4 3; 4 9 16 11];
>> x = P P (B);
>> disp(x)
102 Applied Linear Algebra and Optimization using MATLAB
Program 1.7
MATLAB m-file for Gaussian Elimination by Partial Pivoting
function x=PP(B)
% B = input(0 input matrix in f orm[A/b]0 );
[n, t] = size(B); U = B;
for M = 1:n-1
mx(M ) = abs(U (M, M )); r = M ;
for i = M+1:n
if mx(M ) < abs(U (i, M ))
mx(M)=abs(U(i,M)); r = i; end; end
rw1(1,1:t)=U(r,1:t); rw2(1,1:t)=U(M,1:t);
U(M,1:t)=rw1 ; U(r,1:t)=rw2 ;
for k=M+1:n
m=U(k,M)/U(M,M);
for j=M:t
U (k, j) = U (k, j) − m ∗ U (M, j); end;end
i=n; x(i)=U(i,t)/U(i,i);
for i=n-1:-1:1; s=0;
for k=n:-1:i+1
s = s + U (i, k) ∗ x(k); end
x(i)=(U(i,t)-s)/U(i,i); end; B; U; x; end
1. Suppose we are about to work on the ith column of the matrix. Then
we search that portion of the ith column below and including the di-
agonal and find the element that has the largest absolute value. Let
p denote the index of the row that contains this element.
Total Pivoting
In the case of total pivoting (or complete pivoting), we search for the largest
number (in absolute value) in the entire array instead of just in the first
column, and this number is the pivot. This means that we shall probably
need to interchange the columns as well as rows. When solving a system
of equations using complete pivoting, each row interchange is equivalent to
interchanging two equations, while each column interchange is equivalent
to interchanging the two unknowns.
At the kth step, interchange both the rows and columns of the matrix
so that the largest number in the remaining matrix is used as the pivot
i.e., after the pivoting
|akk | = max|aij |, for i = k, k + 1, . . . , n, j = k, k + 1, . . . , n.
There are times when the partial pivoting procedure is inadequate.
When some rows have coefficients that are very large in comparison to
those in other rows, partial pivoting may not give a correct solution.
Then eliminate the third variable x3 from the second and third rows by
4 1
subtracting the multiples m21 = 16 and m31 = 16 of row 1 from rows 2 and
3, which respectively, gives
For the second elimination step, 1 is the largest absolute coefficient of the
first variable x1 in the second row and third column, so the second and third
columns are interchanged, giving us
Eliminate the first variable x1 from the third row by subtracting the multiple
m32 = 43 of row 2 from row 3, which gives
x1 = 1, x2 = −1, x3 = 1,
Program 1.8
MATLAB m-file for the Gaussian Elimination by Total Pivoting
function x=TP(B)
% B = input(‘input matrix in f orm[A/b]0 );
[n,m]=size(B);U=B; w=zeros(n,n);
for i=1:n; N(i)=i; end
for M = 1:n-1; r=M; c=M;
for i = M:n; for j = M:n
if max(M ) < abs(U (i, j)); max(M)=abs(U(i,j));
r = i; c = j; end; end; end
rw1(1,1:m)=U(r,1:m); rw2(1,1:m)=U(M,1:m);
U(M,1:m)=rw1;U(r,1:m)=rw2 ; cl1(1:n,1)= U(1:n,c);
cl2(1 : n, 1) = U (1 : n, M ); U (1 : n, M ) = cl1(1 : n, 1);
U (1 : n, c) = cl2(1 : n, 1); p = N (M ); N (M ) = N (c);
N (c) = p; w(M, 1 : n) = N ;
for k = M + 1 : n; e = U (k, M )/U (M, M );
for j = M : m; U (k, j) = U (k, j) − e ∗ U (M, j); end; end
i = n; x(i, 1) = U (i, m)/U (i, i);
for i = n − 1 : −1 : 1; s = 0;
for k = n : −1 : i + 1; s = s + U (i, k) ∗ x(k, 1);end
x(i, 1) = (U (i, m) − s)/U (i, i); end
for i=1:n; X(N (i), 1) = x(i, 1); end; B; U ; X; end
>> B = [1 1 1 1; 2 3 4 3; 4 9 16 11];
>> x = T P (B);
>> disp(x)
Total pivoting offers little advantage over partial pivoting and it is signifi-
cantly slower, requiring N = n(n+1)(2n+1)
6
elements to be examined in total.
It is rarely used in practice because interchanging columns changes the
order of the xs and, consequently, add significant and usually unjustified
106 Applied Linear Algebra and Optimization using MATLAB
[A|b] → [I|c],
where I is the identity matrix and c is the column matrix, which represents
the possible solution of the given linear system.
Example 1.34 Solve the following linear system using the Gauss–Jordan
method:
x1 + 2x2 = 3
−x1 − 2x3 = −5
−3x1 − 5x2 + x3 = −4.
Solution. Write the given system in the augmented matrix form
..
1 2 0 . 3
..
.
−1 0 −2 . −5
..
−3 −5 1 . −4
The first elimination step is to eliminate elements a21 = −1 and a31 = −3
by subtracting the multiples m21 = −1 and m31 = −3 of row 1 from rows
Matrices and Linear Systems 107
Program 1.9
MATLAB m-file for the Gauss–Jordan Method
function sol=GaussJ(Ab)
[m,n]=size(Ab);
for i=1:m
Ab(i, :) = Ab(i, :)/Ab(i, i);
for j=1:m
if j == i; continue; end
Ab(j, :) = Ab(j, :) − Ab(j, i) ∗ Ab(i, :);
end; end; sol=Ab;
3. Use the nth row to reduce the nth column to an equivalent identity
matrix column.
4. Repeat step 3 for n–1 through 2 to get the augmented matrix of the
form [I|c].
Matrices and Linear Systems 109
Obviously, the original augmented matrix [A|I] has been transformed to the
augmented matrix of the form [I|A−1 ]. Hence, the solution of the linear
system can be obtained by the matrix multiplication (1.32) as
x1 −0.36 −0.16 0.28 1 1
x2 = 1.6 0.6 −0.8 2 = −2 .
x3 −0.6 −0.2 0.4 6 1.4
A = LU, (1.33)
LU x = b
or can be written as
Ly = b,
where
y = U x.
The unknown elements of matrix L and matrix U are computed by equating
corresponding elements in matrices A and LU in a systematic way. Once
the matrices L and U have been constructed, the solution of system (1.35)
can be computed in the following two steps:
Now defining the three elementary matrices (each of them can be ob-
tained by adding a multiple of row i to row j) associated with these row
operations:
1 0 0 1 0 0 1 0 0
E1 = −2 1 0 , E2 = 0 1 0 , E3 = 0 1 0 .
0 0 1 1 0 1 0 −2 1
Then
1 0 0 1 0 0 1 0 0 1 0 0
E3E2E1 = 0 1 0 0 1 0 −2 1 0 = −2 1 0
0 −2 1 1 0 1 0 0 1 5 −2 1
and
1 0 0 2 4 2 2 4 2
E3E2E1A = −2 1 0 4 9 7 = 0 1 3 = U.
5 −2 1 −2 −2 5 0 0 1
So
A = E1−1 E2−1 E3−1 = LU,
where
1 0 0 1 0 0 1 0 0 1 0 0
E1−1 E2−1 E3−1 = 2 1 0 0 1 0 0 1 0 = 2 1 0 = L.
0 0 1 −1 0 1 0 2 1 −1 2 1
Thus, A = LU is a product of a lower-triangular matrix L and an upper-
triangular matrix U . Naturally, this is called an LU decomposition of A.
Matrices and Linear Systems 115
Doolittle’s Method
In Doolittle’s method (called Gauss factorization), the upper-triangular
matrix U is obtained by forward elimination of the Gaussian elimination
method and the lower-triangular matrix L containing the multiples used in
the Gaussian elimination process as the elements below the diagonal with
unity elements on the main diagonal.
where the unknown elements of matrix L are the used multiples and the
matrix U is the same as we obtained in the forward elimination process.
Example 1.36 Construct the LU decomposition of the following matrix A
by using Gauss factorization (i.e., LU decomposition by Doolittle’s method).
Find the value(s) of α for which the following matrix is
1 −1 α
A = −1 2 −α
α 1 1
singular. Also, find the unique solution of the linear system Ax = [1, 1, 2]T
by using the smallest positive integer value of α.
116 Applied Linear Algebra and Optimization using MATLAB
now we will use only the forward elimination step of the simple Gaussian
elimination method to convert the given matrix A into the upper-triangular
matrix U . Since a11 = 1 6= 0, we wish to eliminate the elements a21 = −1
and a31 = α by subtracting from the second and third rows the appropriate
multiples of the first row. In this case, the multiples are given,
−1 α
m21 = = −1 and m31 = = α.
1 1
Hence,
1 −1 α
0 1 0 .
0 1 + α 1 − α2
(1) (1)
Since a22 = 1 6= 0, we eliminate the entry in the a32 = 1 + α position by
subtracting the multiple m32 = 1+α
1
of the second row from the third row to
get
1 −1 α
0 1 0 .
0 0 1 − α2
Obviously, the original set of equations has been transformed to an upper-
triangular form. Thus,
1 −1 α 1 0 0 1 −1 α
−1 2 −α = −1 1 0 0 1 0 ,
α 1 1 α 1+α 1 0 0 1 − α2
y1 = 1 gives y1 = 1,
−y1 + y2 = 1 gives y2 = 2,
2y1 + 3y2 + y3 = 2 gives y3 = −6.
x1 − x2 + 2x3 = 1 gives x1 = −1
x2 = 2 gives x2 = 2
− 3x3 = −6 gives x3 = 2,
which gives
x1 = −1
x2 = 2
x3 = 2,
the approximate solution of the given system. •
118 Applied Linear Algebra and Optimization using MATLAB
>> A = [1 2 0; −1 0 − 2; −3 − 5 1];
>> B = lu − gauss(A);
>> L = eye(size(B)) + tril(B, −1);
>> U = triu(A);
>> b = [3 − 5 − 4]0 ;
>> y = L \ b;
>> x = U \ y;
Program 1.10
MATLAB m-file for the LU Decomposition Method
function A = lu − gauss(A)
% LU factorization without pivoting
[n,n] = size(A); for i=1:n-1; pivot = A(i,i);
for k=i+1:n; A(k,i)=A(k,i)/pivot;
for j=i+1:n; A(k, j) = A(k, j) − A(k, i) ∗ A(i, j);
end;end; end
There is another way to find the values of the unknown elements of the
matrices L and U , which we describe in the following example.
Solution. Since
1 0 0 u11 u12 u13
A = LU = l21 1 0 0 u22 u23 ,
l31 l32 1 0 0 u33
1 = u11 , u11 = 1
1 = l21 u11 , l21 = 1
2 = l31 u11 , l31 = 2.
2 = u12 , u12 = 2
3 = l21 u12 + u22 , u22 = 3 − 2 = 1
2 = l31 u12 + l32 u22 , l32 = 2 − 4 = −2.
4 = u13 , u13 = 4
3 = l21 u13 + u23 , u23 = 3 − 4 = −1
2 = l31 u13 + l32 u23 + u33 , u33 = 2 − 10 = −8.
Thus, we obtain
1 2 4 1 0 0 1 2 4
1 3 3 = 1 1 0 0 1 1 ,
2 2 2 2 −2 1 0 0 −8
i−1
X
uij = aij − lik ukj , 2 ≤i≤j
k=1
" j−1
#
1 X
lij = aij − lik ukj , i >j≥2
uii . (1.38)
k=1
uij = a1j , i =1
ai1 ai1
lij = = , j =1
u11 a11
1 2 4 x1 −2
0 1 1 x2 = 5 .
0 0 −8 x3 8
which gives
x1 = −6
x2 = 4
x3 = −1,
We can also write the MATLAB m-file called Doolittle.m to get the
solution of the linear system by LU decomposition by using Doolittle’s
method. In order to reproduce the above results using MATLAB com-
mands, we do the following:
>> A = [1 2 4; 1 3 3; 2 2 2];
>> b = [−2 3 − 6];
>> sol = Doolittle(A, b);
122 Applied Linear Algebra and Optimization using MATLAB
Program 1.11
MATLAB m-file for using Doolittle’s Method
function sol = Doolittle(A,b)
[n,n]=size(A); u=A;l=zeros(n,n);
for i=1:n-1; if abs(u(i,i))> 0
for i1=i+1:n; m(i1,i)=u(i1,i)/u(i,i);
for j=1:n
u(i1, j) = u(i1, j) − m(i1, i) ∗ u(i, j);end;end;end;end
for i=1:n; l(i,1)=A(i,1)/u(1,1); end
for j=2:n; for i=2:n; s=0;
for k=1:j-1; s = s + l(i, k) ∗ u(k, j); end
l(i,j)=(A(i,j)-s)/u(j,j); end; end y(1)=b(1)/l(1,1);
for k=2:n; sum=b(k);
for i=1:k-1; sum = sum − l(k, i) ∗ y(i); end
y(k)=sum/l(k,k); end
x(n)=y(n)/u(n,n);
for k=n-1:-1:1; sum=y(k);
for i=k+1:n; sum = sum − u(k, i) ∗ x(i); end
x(k)=sum/u(k,k); end; l; u; y; x
Let D denote the diagonal matrix having the same diagonal elements as
the upper-triangular matrix U ; in other words, D contains the pivots on
its diagonal and zeros everywhere else. Let V be the redefining upper-
triangular matrix obtained from the original upper-triangular matrix U by
dividing each row by its pivot, so that V has all 1s on the diagonal. It
is easily seen that U = DV , which allows any LU decomposition to be
written as
A = LDV,
where L and V are lower- and upper-triangular matrices with 1s on both
of their diagonals. This is called the LDV factorization of A.
L=VT and V = LT .
Note that not every symmetric matrix has an LDLT factorization. How-
ever, if A = LDLT , then A must be symmetric because
Note that
1 3 2 1 0 0
V = 0 1 1 = LT = 3 1 0 .
0 0 1 2 1 1
Thus, we obtain
1 0 0 1 0 0 1 3 2
A = LDLT = 3 1 0 0 −5 0 0 1 1 ,
2 1 1 0 0 3 0 0 1
Crout’s Method
Crout’s method, in which matrix U has unity on the main diagonal, is
similar to Doolittle’s method in all other aspects. The L and U matrices
are obtained by expanding the matrix equation A = LU term by term to
determine the elements of the L and U matrices.
Solution. Since
l11 0 0 1 u12 u13
A = LU = l21 l22 0 0 1 u23 ,
l31 l32 l33 0 0 1
performing the multiplication on the right-hand side gives
1 2 3 l11 l11 u12 l11 u13
6 5 4 = l21 l21 u12 + l22 l21 u13 + l22 u23 .
2 5 6 l31 l31 u12 + l32 l31 u13 + l32 u23 + l33
Then equate elements of the first column to obtain
1 = l11
6 = l21
2 = l31 .
Then equate elements of the second column to obtain
2 = l11 u12 , u12 = 2
j−1
X
lij = aij − lik ukj , i ≥ j, i = 1, 2, . . . , n
k=1
i−1
1 X
uij = [aij − lik ukj ], i < j, j = 2, 3, . . . , n
lii . (1.39)
k=1
lij = ai1 , j=1
aij
uij = , i=1
a11
y1 = 1 gives y1 = 1
6y1 − 7y2 = −1 gives y2 = 1
2y1 + y2 − 2y3 = 5 gives y3 = −1.
1 2 3 x1 1
0 1 2 x2 = 1 .
0 0 1 x3 −1
>> A = [1 2 3; 6 5 4; 2 5 6];
>> b = [1 − 1 5];
>> sol = Crout(A, b);
Matrices and Linear Systems 129
Program 1.12
MATLAB m-file for the Crout’s Method
function sol = Crout(A, b)
[n,n]=size(A); u=zeros(n,n); l=u;
for i=1:n; u(i,i)=1; end
l(1,1)=A(1,1);
for i=2:n
u(1,i)=A(1,i)/l(1,1); l(i,1)=A(i,1); end
for i=2:n; for j=2:n; s=0;
if i <= j; K=i-1;
else; K=j-1; end
for k=1:K; s = s + l(i, k) ∗ u(k, j); end
if j > i; u(i,j)=(A(i,j)-s)/l(i,i); else
l(i,j)=A(i,j)-s; end;end;end
y(1)=b(1)/l(1,1);
for k=2:n; sum=b(k);
for i=1:k-1; sum = sum − l(k, i) ∗ y(i); end
y(k)=sum/l(k,k); end
x(n)=y(n)/u(n,n);
for k=n-1:-1:1; sum=y(k);
for i=k+1:n; sum = sum − u(k, i) ∗ x(i); end
x(k)=sum/u(k,k); end; l; u; y; x;
Example 1.43 Find the determinant and inverse of the following matrix
using LU decomposition by Doolittle’s method:
1 −2 1
A = 1 −1 1 .
1 1 2
Solution. We know that
1 −2 1 1 0 0 u11 u12 u13
A = 1 −1 1 = m21 1 0 0 u22 u23 = LU.
1 1 2 m31 m32 1 0 0 u33
Matrices and Linear Systems 131
Now we will use only the forward elimination step of the simple Gaussian
elimination method to convert the given matrix A into the upper-triangular
matrix U . Since a11 = 1 6= 0, we wish to eliminate the elements a21 = 1
and a31 = 1 by subtracting from the second and third rows the appropriate
multiples of the first row. In this case, the multiples are given as
Hence,
1 −2 1
0 1 0 .
0 3 1
(1) (1)
Since a22 = 1 6= 0, we eliminate the entry in the a32 = 3 position by
subtracting the multiple m32 = 3 of the second row from the third row to
get
1 −2 1
0 1 0 .
0 0 1
Obviously, the original set of equations has been transformed to an upper-
triangular form. Thus,
1 −2 1 1 0 0 1 −2 1
1 −1 1 = 1 1 0 0 1 0 ,
1 1 2 1 3 1 0 0 1
To find the inverse of matrix A, first we will compute the inverse of the
lower-triangular matrix L−1 from
0
1 0 0 l11 0 0 1 0 0
LL−1 = 1 1 0 l21 0 0
l22 0 = 0 1 0 =I
0 0 0
1 3 1 l31 l32 l33 0 0 1
132 Applied Linear Algebra and Optimization using MATLAB
can be obtained
0 0
l22 = 1, l32 = −3.
Finally, the solution of the third linear system
1 0 0 0 0
1 1 0 0 = 0
0
1 3 1 l33 1
0
gives l33 = 1.
Hence, the elements of the matrix L−1 are
1 0 0
L−1 = −1 1 0 ,
2 −3 1
For LU decomposition we have not used pivoting for the sake of sim-
plicity. However, pivoting is important for the same reason as in Gaussian
elimination. We know that pivoting in Gaussian elimination is equivalent
to interchanging the rows of the coefficients matrix together with the terms
on the right-hand side. This indicates that pivoting may be applied to LU
decomposition as long as the interchanging is applied to the left and right
terms in the same way. When performing pivoting in LU decomposition,
the changes in the order of the rows are recorded. The same reordering is
then applied to the right-hand side terms before starting the solution in
accordance with the forward elimination and backward substitution steps.
Indirect LU Decomposition
and replacing the matrix A by P A. For example, using the above matrix
A, we have
1 0 0 2 2 −4 2 2 −4
P A = 0 0 1 2 2 −1 = 3 2 −3 .
0 1 0 3 2 −3 2 2 −1
From this multiplication we see that rows 2 and 3 of the original matrix A
are interchanged, and the resulting matrix P A has a LU factorization and
we have
1 0 0 2 2 −4 2 2 −4
1.5 1 0 0 −1 3 = 3 2 −3 .
1 0 1 0 0 3 2 2 −1
A0 = P A, equivalently A = P −1 A0 ,
A0 = P A = LU
and so
A = P −1 LU = (P T L)U
since P −1 = P T . The determinant of A may now be written as
or
det(A) = β det(L) det(U ),
where β = det(P −1 ) equals −1 or +1 depending on whether the number
pivoting is odd or even, respectively. •
One can use the MATLAB built-in lu function to obtain the permuta-
tion matrix P so that the P A matrix has a LU decomposition:
>> A = [0 1 2; −1 4 2; 2 2 1];
>> [L, U, P ] = lu(A);
It will give us the permutation matrix P and the matrices L and U as
follows:
0 0 1
P = 0 1 0
1 0 0
and
1 0 0 2 2 1
P A = −0.5 1 0 0 5 2.5 = LU.
0 0.2 1 0 0 1.5
So
A = P −1 LU
Matrices and Linear Systems 137
or
0 0.2 1 2 2 1
A = (P T L)U = −0.5 1 0 0 5 2.5 .
1 0 0 0 0 1.5
then:
1. Show that A does not have LU factorization;
Note that during this elimination process two row interchanges were needed,
which means we got two elementary permutation matrices of the inter-
changes (from Theorem 1.23), which are
0 0 1 1 0 0
p1 = 0 1 0 and p2 = 0 0 1 .
1 0 0 0 1 0
3
m21 = 2, m31 = 0, m32 = −
2
as follows:
3 2 5 3 2 5 3 2 5
6 2 4 → 6 −2 −6 → 0 −2 −6 = U.
0 3 3 0 3 3 0 0 −6
Matrices and Linear Systems 139
Thus, P A = LU , where
1 0 0 3 2 5
L= 2 1 0 and U = 0 −2 −6 .
0 −3/2 1 0 0 −6
(3) Solve the first system Ly = P b = [4, 3, 6]T for unknown vector y, i.e.,
1 0 0 y1 4
2 1 0 y2 = 3 .
0 −3/2 1 y3 6
y1 = 4 gives y1 = 4
2y1 + y2 = 3 gives y2 = −5
− 3/2y2 + y3 = 6 gives y3 = −1.5.
Then solve the second system U x = y for the unknown vector x, i.e.,
3 2 5 x1 4
0 −2 −6 x2 = −5 .
0 0 −6 x3 −1.5
The function
a11 a12 · · · a1n
x 1
a21 a22 · · · a2n x2
xT Ax = x1 x2 · · · xn
a31 a32 · · · a3n
.
.. .. .. .. ..
. . . .
xn
an1 an2 · · · ann
or n X
n
X
xT Ax = aij xi xj (1.42)
i=1 j=1
can be used to represent any quadratic polynomial in the variables x1 , x2 , . . . , xn
and is called a quadratic form. A matrix is said to be positive-definite if
its quadratic form is positive for all real nonzero vectors x, i.e.,
xT Ax > 0, for every n-dimensional column vector x 6= 0.
Example 1.45 The matrix
4 −1 0
A = −1 4 −1
0 −1 4
is positive-definite and suppose x is any nonzero three-dimensional column
vector, then
4 −1 0 x1
xT Ax = x1 x2 x3 −1
4 −1 x2
0 −1 4 x3
Matrices and Linear Systems 141
or
4x1 − x2
= x1 x2 x3 −x1 + 4x2 − x3 .
− x2 + 4x3
Thus,
xT Ax = 4x21 − 2x1 x2 + 4x22 − 2x2 x3 + 4x23 .
After rearranging the terms, we have
Hence,
3x21 + (x1 − x2 )2 + 2x22 + (x2 − x3 )2 + 3x23 > 0,
unless x1 = x2 = x3 = 0. •
The principal minors of a matrix A are the square submatrices lying in the
upper-left hand corner of A. An n × n matrix A has n of these principal
minors. For example, for the matrix
6 2 1
A = 2 3 2 ,
1 1 2
142 Applied Linear Algebra and Optimization using MATLAB
det(6) = 6 > 0,
6 2
det = 18 − 4 = 14 > 0,
2 3
6 2 1
det 2 3 2 = 19 > 0.
1 1 2
det(4) = 4 > 0,
4 −1
det = 12 − 1 = 11 > 0,
−1 3
4 −1 2
det −1 3 0 = 43 > 0.
2 0 5
we can have
3 4 4
AT A = 4 6 5 .
4 5 6
det(3) = 3 > 0,
3 4
det = 18 − 16 = 2 > 0,
4 6
3 4 4
det 4 6 5 = 1 > 0.
4 5 6
Cholesky Method
The Cholesky method (or square root method) is of the same form as
Doolittle’s method and Crout’s method except it is limited to equations
involving symmetrical coefficient matrices. In the case of a symmetric
and positive-definite matrix A it is possible to construct an alternative
triangular factorization with a saved number of calculations compared with
previous factorizations. Here, we decompose the matrix A into the product
of LLT , i.e.,
A = LLT , (1.43)
144 Applied Linear Algebra and Optimization using MATLAB
2. Solve LT x = y, for x.
(using backward substitution)
In this procedure, it is necessary to take the square root of the elements
on the main diagonal of the coefficient matrix. However, for a positive-
definite matrix the terms on its main diagonal are positive, so no difficulty
will arise when taking the square root of these terms.
Example 1.46 Construct the LU decomposition of the following matrix
using the Cholesky method:
1 1 2
A = 1 2 4 .
2 4 9
Solution. Since
l11 0 0 l11 l21 l31
A = LLT = l21 l22 0 0 l22 l32 ,
l31 l32 l33 0 0 l33
2
√
1 = l11 gives l11 = 1=1
√
Note that l11 could be − 1 and so the matrix L is not (quite) unique.
Now equate elements of the second column to obtain
2 2
2 = l21 + l22 gives l22 = 1
4 = l31 l21 + l32 l22 gives l32 = 2.
2 2 2
9 = l31 + l32 + l33 gives l33 = 1.
Thus, we obtain
1 1 2 1 0 0 1 1 2
1 2 4 = 1 1 0 0 1 2 ,
2 4 9 2 2 1 0 0 1
The method fails if ljj = 0 and the expression inside the square root is
negative, in which case all of the elements in column j are purely imaginary.
There is, however, a special class of matrices for which these problems don’t
occur.
The Cholesky method provides a convenient method for investigat-
ing the positive-definiteness of symmetric matrices. The formal definition
xT Ax > 0, for all x 6= 0, is not easy to verify in practice. However, it is
relatively straightforward to attempt the construct of a Cholesky decom-
position of a symmetric matrix.
Solution. Since
l11 0 0 l11 l21 l31
A = LLT = l21 l22 0 0 l22 l32 ,
l31 l32 l33 0 0 l33
performing the multiplication on the right-hand side gives
2
9 3 6 l11 l11 l21 l11 l31
3 10 8 = l11 l21 l21 2 2
+ l22 l21 l31 + l22 l32 .
2 2 2
6 8 9 l11 l31 l31 l21 + l22 l32 l31 + l32 + l33
Then equate elements of the first column to obtain
2
√
9 = l11 gives l11 = 9=3
y1 = 2, y1 = 2
y1 + y2 = 1, y2 = −1
2y1 + 2y2 + y3 = 1, y3 = −1.
x1 + x2 + 2x3 = 2 gives x1 = 3
x2 + 2x3 = −1 gives x2 = 1
x3 = −1 gives x3 = −1,
Now use the following MATLAB commands to obtain the above results:
>> A = [1 1 2; 1 2 4; 2 4 9];
>> b = [2 1 1];
>> sol = Cholesky(A, b);
Example 1.49 Find the bounds on α for which the Cholesky factorization
of the following matrix with real elements
1 2 α
A = 2 8 2α
α 2α 9
is possible.
Solution. Since
l11 0 0 l11 l21 l31
A = LLT = l21 l22 0 0 l22 l32 ,
l31 l32 l33 0 0 l33
150 Applied Linear Algebra and Optimization using MATLAB
2
1 2 α l11 l11 l21 l11 l31
2
2 8 2α = l11 l21 l21 2
+ l22 l21 l31 + l22 l32 .
2 2 2
α 2α 9 l11 l31 l31 l21 + l22 l32 l31 + l32 + l33
2
√
1 = l11 gives l11 = 1=1
√
Note that l11 could be − 1 and so matrix L is not (quite) unique.
Now equate elements of the second column to obtain
2 2
8 = l21 + l22 gives l22 = 2
2α = l31 l21 + l32 l22 gives l32 = 0.
2 2 2
√
9 = l31 + l32 + l33 gives l33 = 9 − α2 ,
Program 1.13
MATLAB m-file for the Cholesky Method
function sol = Cholesky(A, b)
[n,n]=size(A); l=zeros(n,n); u=l;
l(1, 1) = (A(1, 1)) \ 0.5; u(1,1)=l(1,1);
for i=2:n; u(1,i)=A(1,i)/l(1,1);
l(i,1)=A(i,1)/u(1,1); end
for i=2:n; for j=2:n; s=0;
if i <= j; K=i-1; else; K=j-1; end
for k=1:K; s = s + l(i, k) ∗ u(k, j); end
if j > i; u(i,j)=(A(i,j)-s)/l(i,i);
elseif i == j
l(i, j) = (A(i, j) − s) \ 0.5; u(i,j)=l(i,j);
else; l(i,j)=(A(i,j)-s)/u(j,j); end; end; end
y(1)=b(1)/l(1,1);
for k=2:n; sum=b(k);
for i=1:k-1; sum = sum − l(k, i) ∗ y(i); end
y(k)=sum/l(k,k); end
x(n)=y(n)/u(n,n);
for k=n-1:-1:1; sum=y(k);
for i=k+1:n; sum = sum − u(k, i) ∗ x(i); end
x(k)=sum/u(k,k); end; l; u; y; x;
Solution. By using the simple Gauss elimination method, one can convert
the given matrix into the upper-triangular matrix
4 −2 4
U = 0 1 4
0 0 9
152 Applied Linear Algebra and Optimization using MATLAB
where
1 0 0 2 0 0 2 0 0
L̂ = LD1/2 = −0.5 1 0 0 1 0 = −1 1 0
1 4 1 0 0 3 2 4 3
Matrices and Linear Systems 153
and
2 0 0 1 −0.5 1 2 −1 2
L̂T = D1/2 V = 0 1 0 0 1 4 = 0 1 4 .
0 0 3 0 0 1 0 0 3
Thus, we obtain
2 0 0 2 −1 2
A = L̂L̂T = −1 1 0 0 1 4 ,
2 4 3 0 0 3
1. Matrix A is nonsingular.
Example 1.53 Solve the following linear system using the simple Gaus-
sian elimination method and also find the LU decomposition of the matrix
using Doolittle’s method and Crout’s method:
5x1 + x2 + x3 = 7
2x1 + 6x2 + x3 = 9
x1 + 2x2 + 9x3 = 12.
and since a11 = 5 6= 0, we can eliminate the elements a21 and a31 by
subtracting from the second and third rows the appropriate multiples of the
first row. In this case the multiples are given,
2 1
m21 = = 0.4 and m31 = = 0.2.
5 5
Hence,
..
5 1 1 . 7
0 5.6 0.6 ... 6.2 .
..
0 1.8 8.8 . 10.6
(1) (1)
Since a22 = 5.6 6= 0, we eliminate the entry in the a32 position by sub-
1.8
tracting the multiple m32 = 5.6 = 0.32 of the second row from the third row
156 Applied Linear Algebra and Optimization using MATLAB
to get
.
5 1 ..
1 7
0 5.6 0.6 ... 6.2 .
..
0 0 8.6 . 8.6
Obviously, the original set of equations has been transformed to an upper-
triangular form. All the diagonal elements of the obtaining upper-triangular
matrix are nonzero, which means that the coefficient matrix of the given
system is nonsingular, therefore, the given system has a unique solution.
Now expressing the set in algebraic form yields
5x1 + x2 + x3 = 7
5.6x2 + 0.6x3 = 6.2
8.6x3 = 8.6.
Now use backward substitution to get the solution of the system as
8.6x3 = 8.6 gives x3 = 1
5.6x2 = −0.6x3 + 6.2 = −0.6 + 6.2 = 5.6 gives x2 = 1
5x1 = 7 − x2 − x3 = 7 − 1 − 1 = 5 gives x1 = 1.
We know that when using LU decomposition by Doolittle’s method the un-
known elements of matrix L are the multiples used and the matrix U is the
same as we obtained in the forward elimination process of the simple Gauss
elimination. Thus, the LU decomposition of matrix A can be obtained by
using Doolittle’s method as follows:
5 1 1 1 0 0 5 1 1
A = 2 6 1 0.4 1 0 0 5.6 0.6 = LU.
1 2 9 0.5 0.32 1 0 0 8.6
Similarly, the LU decomposition of matrix A can be obtained by using
Crout’s method as
5 1 1 5 0 0 1 0.2 0.2
A= 2 6 1 2 5.6 0 0 1 0.1 = LU.
1 2 9 1 1.8 8.6 0 0 1
Thus, the conditions of Theorem 1.30 are satisfied. •
Matrices and Linear Systems 157
After finding the values for li and ui , then they are used along with the
elements ci , to solve the tridiagonal system (1.47) by solving the first bidi-
agonal system
Ly = b, (1.50)
for y by using forward substitution,
y 1 = b1
, (1.51)
yi = bi − li yi−1 , i = 2, 3, . . . , n
U x = y, (1.52)
The entire process for solving the original system (1.47) requires 3n ad-
ditions, 3n multiplications, and 2n divisions. Thus, the total number of
multiplications and divisions is approximately 5n.
x1 + x2 = 1
x1 + 2x2 + x3 = 0
x2 + 3x3 + x4 = 1
x3 + 4x4 = 1.
Then the elements of the L and U matrices can be computed by using (1.48)
as follows:
u1 = α1 = 1
β2 1
l2 = = =1
u1 1
u2 = α2 − l2 c1 = 2 − (1)1 = 1
b3 1
l3 = = =1
u2 1
u3 = α3 − l3 c2 = 3 − (1)1 = 2
b4 1
l4 = =
u3 2
1 7
u4 = α4 − l4 c3 = 4 − ( )1 = .
2 2
After finding the elements of the bidiagonal matrices L and U , we solve the
160 Applied Linear Algebra and Optimization using MATLAB
1 0 0 0 y1 1
1 1 0 0 y2 0
= .
0 1 1 0 y3 1
0 0 21 1 y4 1
1 1 0 0 x1 1
0 1 1 0 x2 −1
=
2 .
0 0 2 1 x3
0 0 0 27 x4 0
Program 1.14
MATLAB m-file for LU Decomposition for a Tridiagonal System
function sol=TRiDLU(Tb)
[m,n]=size(Tb); L=eye(m); U=zeros(m);
U(1,1)=Tb(1,1);
for i=2:m
U (i − 1, i) = T b(i − 1, i);
L(i, i − 1) = T b(i, i − 1)/U (i − 1, i − 1);
U (i, i) − L(i, i − 1) ∗ T b(i − 1, i); end
disp(’The lower-triangular matrix’) L;
disp(’The upper-triangular matrix’) U;
y = inv(L) ∗ T b(:, n); x = inv(U ) ∗ y;
Vector Norms
There are three norms in Rn that are most commonly used in applications,
called l1 -norm, l2 -norm, and l∞ -norm, and are defined for the given vectors
Matrices and Linear Systems 163
x = [x1 , x2 , . . . , xn ]T as
n
X
kxk1 = |xi |
i=1
n
!1/2
X
kxk2 = x2i
i=1
The l1 -norm is called the absolute norm, the l2 -norm is frequently called
the Euclidean norm as it is just the formula for distance in ordinary three-
dimensional Euclidean space extended to dimension n, and finally, the
l∞ -norm is called the maximum norm or occasionally the uniform norm.
All these three norms are also called the natural norms.
Matrix Norms
A matrix norm is a measure of how well one matrix approximates another,
or, more accurately, of how well their difference approximates the zero ma-
trix. An iterative procedure for inverting a matrix produces a sequence
of approximate inverses. Since, in practice, such a process must be termi-
nated, it is desirable to have some measure of the error of an approximate
inverse.
1. kAk > 0, A 6= 0;
2. kAk = 0, A = 0;
5. kA + Bk ≤ kAk + kBk;
Matrices and Linear Systems 165
6. kABk ≤ kAkkBk;
7. kA − Bk ≥ kAk − kBk.
Several norms for matrices have been defined, and we shall use the following
three natural norms l1 , l2 , and l∞ for a square matrix of order n:
n
!
X
kAk1 = max |aij | = maximum column sum,
j
i=1
n
!
X
kAk∞ = max |aij | = row sum norm.
i
j=1
The l1 -norm and l∞ -norm are widely used because they are easy to cal-
culate. The matrix norm kAk2 that corresponds to the l2 -norm is related
to the eigenvalues of the matrix. It sometimes has special utility because
no other norm is smaller than this norm. It, therefore, provides the best
measure of the size of a matrix, but is also the most difficult to compute.
We will discuss this natural norm later in the chapter.
m X
n
!1/2
X
kAkF = |aij |2 .
i=1 j=1
so
kAk1 = max{8, 9, 10} = 10.
Also,
3
X
|a1j | = |4| + |2| + | − 1| = 7,
j=1
3
X
|a2j | = |3| + |5| + | − 2| = 10,
j=1
X3
|a3j | = |1| + | − 2| + |7| = 10,
j=1
so
kAk∞ = max{7, 10, 10} = 10.
Finally, we have
>> A = [4 2 − 1; 3 5 − 2; 1 − 2 − 7];
>> B = norm(A, 1)
B=
10
The l∞ -norm of a matrix A is:
>> A = [4 2 − 1; 3 5 − 2; 1 − 2 − 7];
>> B = norm(A, inf )
B=
10
Finally, the Frobenius norm of the matrix A is:
>> A = [4 2 − 1; 3 5 − 2; 1 − 2 − 7];
>> B = norm(A,0 f ro0 )
B=
10.6301
Program 1.15
MATLAB m-file for Finding the Residual Vector
function r=RES(A,b,x0)
[n,n]=size(A);
for i=1:n; R(i) = b(i);
for j=1:n
R(i)=R(i)-A(i,j)*x0(j);end
RES(i)=R(i); end
r=RES’
has the approximate solution x∗ = [3, 0]T . To see how good this solution
is, we compute the residual, r = [0, −0.0002]T .
has the exact solution x = [1, 1, 1, 1]T and the approximate solution due to
Gaussian elimination without pivoting is
We found that all the elements of the residual for the second case (with
pivoting) are less than 0.6 × 10−7 , whereas for the first case (without piv-
oting) they are as large as 0.2 × 10−4 . Even without knowing the exact
solution, it is clear that the solution obtained in the second case is much
better than the first case. The residual provides a reasonable measure of
the accuracy of a solution in those cases where the error is primarily due
to the accumulation of round-off errors.
We have seen that for ill-conditioned systems the residual is not neces-
sarily a good measure of the accuracy of a solution. How then can we tell
when a system is ill-conditioned? In the following we discuss some possible
indicators of ill-conditioned systems.
Note that the condition number K(A) for A depends on the matrix norm
used and can, for some matrices, vary considerably as the matrix norm is
changed. Since
The condition numbers provide bounds for the sensitivity of the solution
of a set of equations to changes in the coefficient matrix. Unfortunately,
the evaluation of any of the condition numbers of a matrix A is not a trivial
task since it is necessary first to obtain its inverse.
Example 1.57 Compute the condition number of the following matrix us-
ing the l∞ -norm:
2 −1 0
A = 2 −4 −1 .
−1 0 2
Solution. The condition number of a matrix is defined as
and
n 8 −2 −1 3 −4 −2 4 −1 6 o
kA−1 k∞ = max + + , + + , + + ,
13 13 13 13 13 13 13 13 13
which gives
11
kA−1 k∞ = .
13
172 Applied Linear Algebra and Optimization using MATLAB
Therefore,
−1 11
K(A) = kAk∞ kA k∞ = (7) ≈ 5.9231.
13
Depending on the application, we might consider this number to be rea-
sonably small and conclude that the given matrix A is reasonably well-
conditioned. •
To get the above results using MATLAB commands, we do the follow-
ing:
>> A = [2 − 1 0; 2 − 4 − 1; −1 0 2];
>> Ainv = inv(A)
>> K(A) = norm(A, inf ) ∗ norm(Ainv, inf );
K(A) =
5.9231
Some matrices are notoriously ill-conditioned. For example, consider the
4 × 4 Hilbert matrix
1 1 1
1
2 3 4
1 1 1 1
2 3 4 5
H= ,
1 1 1 1
3 4 5 6
1 1 1 1
4 5 6 7
whose entries are defined by
1
aij = , for i, j = 1, 2, . . . , n.
(i + j − 1)
The inverse of the matrix H can be obtained as
16 −120 240 −140
−120 1200 −2700 1680
H −1 =
240 −2700
.
6480 −4200
−140 1680 −4200 2800
Matrices and Linear Systems 173
which is quite large. Note that the condition numbers of Hilbert matri-
ces increase rapidly as the sizes of the matrices increase. Therefore, large
Hilbert matrices are considered to be extremely ill-conditioned.
for which det A = 10−14 ≈ 0. One can easily find the condition number of
the given matrix as
The condition number of a matrix K(A) using the l2 -norm can be com-
puted by the built-in function cond command in MATLAB as follows:
>> A = [1 − 1 2; 3 1 − 1; 2 0 1];
>> K(A) = cond(A);
K(A) =
19.7982
kx − x∗ k ≤ krkkA−1 k, (1.56)
174 Applied Linear Algebra and Optimization using MATLAB
Ax − Ax∗ = b − (b − r) = r,
kbk
kbk ≤ kAkkxk, or kxk ≥ .
kAk
Hence,
kx − x∗ k kA−1 kkrk krk
≤ ≤ K(A) .
kxk kbk/kAk kbk
The inequalities (1.56) and (1.57) imply that the quantities kA−1 k and
K(A) can be used to give an indication of the connection between the
residual vector and the accuracy of the approximation. If the quantity
K(A) ≈ 1, the relative error will be fairly close to the relative residual.
But if K(A) >> 1, then the relative error could be many times larger than
the relative residual.
x1 + x2 − x3 = 1
x1 + 2x2 − 2x3 = 0
−2x1 + x2 + x3 = −1.
Matrices and Linear Systems 175
>> A = [1 1 − 1; 1 2 − 2; −2 1 1];
>> K(A) = norm(A, inf ) ∗ norm(inv(A), inf );
(b) The residual vector can be calculated as
r = b − Ax∗
1 1 1 −1 2.01
= 0 − 1 2 −2 1.01 .
−1 −2 1 1 1.98
176 Applied Linear Algebra and Optimization using MATLAB
and it gives
krk∞ = 0.07.
>> A = [1 1 − 1; 1 2 − 2; −2 1 1];
>> b = [1 0 − 1]0 ;
>> x0 = [2.01 1.01 1.98]0 ;
>> r = RES(A, b, x0);
>> rnorm = norm(r, inf );
(c) From (1.57), we have
kx − x∗ k krk
≤ K(A) .
kxk kbk
By using parts (a) and (b) and the value kbk∞ = 1, we obtain
kx − x∗ k (0.07)
≤ (22.5) = 1.575.
kxk 1
After applying the forward elimination step of the simple Gauss elimination
method, we obtain
..
1 1 −1 . −0.04
0 1 −1 ... −0.03 .
..
0 0 2 . 0.04
Now by using backward substitution, we obtain the solution
Conditioning
Let us consider the conditioning of the linear system
Ax = b. (1.59)
Case 1.1 Suppose that the right-hand side term b is replaced by b + δb,
where δb is an error in b. If x + δx is the solution corresponding to the
right-hand side b + δb, then we have
Ax + Aδx = b + δb,
Aδx = δb.
δx = A−1 δb.
Thus, the change kδxk in the solution is bounded by kA−1 k times the change
kδbk in the right-hand side.
The conditioning of the linear system is connected with the ratio between
kδxk kδbk
the relative error and the relative change in the right-hand side,
kxk kbk
which gives
Thus, the relative change in the solution is bounded by the condition num-
ber of the matrix times the relative change in the right-hand side. When
the product in the right-hand side is small, the relative change in the so-
lution is small.
or
Aδx = −δA(x + δx).
Multiplying by A−1 , we get
or
kδxk(1 − kA−1 kkδAk) ≤ kA−1 kkδAkkxk,
which can be written as
kδxk kA−1 kkδAk K(A)kδAk/kAk
≤ −1
= . (1.64)
kxk (1 − kA kkδAk) (1 − kA−1 kkδAk)
If the product kA−1 kkδAk is much smaller than 1, the denominator in
(1.64) is near 1. Consequently, when kA−1 kkδAk is much smaller than 1,
then (1.64) implies that the relative change in the solution is bounded by
the condition number of a matrix times the relative change in the coefficient
matrix.
Case 1.3 Suppose that there is a change in the coefficient matrix A and the
right-hand side term b together, and if x + δx is the solution corresponding
to the coefficient matrix A + δA and the right-hand side b + δb, then we
have
(A + δA)(x + δx) = (b + δb), (1.65)
which implies that
Ax + Aδx + xδA + δAδx = b + δb
or
Aδx + δxδA = (δb − xδA).
Multiplying by A−1 , we get
δx(I + A−1 δA) = A−1 (δb − xδA)
or
δx = (I + A−1 δA)−1 A−1 (δb − xδA). (1.66)
Since we know that if A is nonsingular and δA is the error in A, we obtain
kA−1 δAk ≤ kA−1 kkδAk < 1, (1.67)
it then follows that (see Fröberg 1969) the matrix (I+A−1 δA) is nonsingular
and
1 1
k(I + A−1 δA)−1 k ≤ −1
≤ −1
. (1.68)
1 − kA δAk 1 − kA kkδAk
180 Applied Linear Algebra and Optimization using MATLAB
1.6 Applications
In this section we discuss applications of linear systems. Here, we will solve
or tackle a variety of real-life problems from several areas of science.
Solution. Observe that in this example we are given three points and we
want to find a polynomial of degree 2 (one less than the number of data
points). Let the polynomial be
p(x) = a0 + a1 x + a2 x2 .
We are given three points and shall use these three sets of information to
determine the three unknowns a0 , a1 , and a2 . Substituting
x = 1, y = 6; x = 2, y = 3; x = 3, y = 2,
in turn, into the polynomial leads to the following system of three linear
equations in a0 , a1 , and a2 :
a0 + a1 + a2 = 6
a0 + 2a1 + 4a2 = 3
a0 + 3a1 + 9a2 = 2.
Solve this system for a2 , a1 , and a0 using the Gauss elimination method:
. . .
1 1 1 .. 6 1 1 1 .. 6 1 1 1 .. 6
1 2 4 ... 0 1 3 ... −3 ≈ 0 1 3 ... −3 .
3
≈
. .. ..
1 3 9 .. 2 0 2 8 . −4 0 0 2 . 2
Now use backward substitution to get the solution of the system (Fig-
ure 1.4),
2a2 = 2 gives a2 = 1
a1 + 3a2 = −3 gives a1 = −6
a0 + a1 + a2 = 6 gives a0 = 11.
Thus,
p(x) = 11 − 6x + x2
is the required the polynomial. •
Matrices and Linear Systems 183
1. Junctions: All the current flowing into a junction must flow out of
it.
2. Paths: The sum of the IR terms (where I denotes current and R
resistance) in any direction around a closed path is equal to the total voltage
in the path in that direction. •
Example 1.60 Consider the electric network in Figure 1.5. Let us deter-
mine the currents through each branch of this network.
Solution. The batteries are 8 volts and 16 volts. The resistances are 1
ohm, 4 ohms, and 2 ohms. The current entering each battery will be the
184 Applied Linear Algebra and Optimization using MATLAB
Junctions
Junction B : I1 + I2 = I3
Junction D : I3 = I1 + I2
These two equations result in a single linear equation
I1 + I2 − I3 = 0.
Paths
The problem thus reduces to solving the following system of three linear
equations in three variables I1 , I2 , and I3 :
I1 + I2 − I3 = 0
4I1 + I3 = 8
4I2 + I3 = 16.
Solve this system for I1 , I2 , and I3 using the Gauss elimination method:
.. .. .
1 1 −1 . 0 1 1 −1 . 0 1 1 −1 .. 0
. . .
4
0 1 .. 8 ≈
0 −4 5 .. ≈ 0 −4
8 5 .. 8 .
. . .
0 4 1 .. 16 0 4 1 .. −4 0 0 6 .. 24
Now use backward substitution to get the solution of the system:
6I3 = 24 gives I3 = 4
I1 + I2 − I3 = 0 gives I1 = 1.
Thus, the currents are I1 = 1, I2 = 3, and I3 = 4. The units are amps.
The solution is unique, as is to be expected in this physical situation. •
Traffic Flow
.
1 0 0 0 0 1 −1 .. 600
.
0 1 0 0 0 0 −1 .. 500
..
0 0 1 0 0 0 −1 . −200
≈ ..
.
1 −1 .
0 0 0 1 0 600
..
0 0 0 0 1 −1 1 . 000
..
0 0 0 0 0 0 0 . 000
The system of equations that corresponds to this form is:
Let us now examine what the flows in the other branches will be when
this minimum flow along Adams is attained, when x7 gives
x1 = −x6 + 800
x2 = + 700
x3 = 000
x4 = −x6 + 800
x5 = x6 − 200.
Since x7 = 200 implies that x3 = 0 and vice-versa, we see that the minimum
flow in branch x7 can be attained by making x3 = 0; i.e., by closing branch
DE to traffic. •
Suppose we have a thin rectangular metal plate whose edges are kept at
fixed temperatures. As an example, let the left edge be 0o C , the right edge
2o C, and the top and bottom edges 1o C (Figure 1.7). We want to know
the temperature inside the plate. There are several ways of approaching
this kind of problem. The simplest approach of interest to us will be
the following type of approximation: we shall overlay our plate with finer
and finer grids, or meshes. The intersections of the mesh lines are called
mesh points. Mesh points are divided into boundary and interior points,
depending on whether they lie on the boundary or the interior of the plate.
We may consider these points as heat elements, such that each influences
its neighboring points. We need the temperature of the interior points,
190 Applied Linear Algebra and Optimization using MATLAB
given the temperature of the boundary points. It is obvious that the finer
the grid, the better the approximation of the temperature distribution of
the plate. To compute the temperature of the interior points, we use the
following principle.
Theorem 1.34 (Mean Value Property for Heat Conduction)
1
x2 = (x1 + x4 + 3)
4
1
x3 = (x1 + x4 + 1)
4
1
x4 = (x2 + x3 + 3) .
4
The problem thus reduces to solving the following system of four linear
equations in four variables x1 , x2 , x3 , and x4 :
4x1 − x2 − x3 = 1
−x1 + 4x2 − x4 = 3
−x1 + 4x3 − x4 = 1
− x2 − x3 + 4x4 = 3.
Solve this system for x1 , x2 , x3 , and x4 using the Gauss elimination method:
..
4 −1 −1 0 . 1
..
4 −1 −1 0 . 1 15 1 . 13
..
. 0 − −1
−1 4 0 −1 .. 3 4 4 4
≈ ··· ≈ .
−1 .
.. 1 56 16 .
. 22
0 4 −1 0 0 − .
15 15 15
..
0 −1 −1 4 . 3 24 .. 30
0 0 0 .
7 7
Now use backward substitution to get the solution of the system:
24 30 5
x4 = gives x4 =
7 7 4
56 16 22 3
x3 − x4 = gives x3 =
15 15 15 4
15 1 13 5
x2 − x3 − x4 = gives x2 =
4 4 4 4
3
4x1 − x2 − x3 = 1 gives x1 = .
4
Thus, the temperatures are x1 = 34 , x2 = 54 , x3 = 34 , and x4 = 54 . •
192 Applied Linear Algebra and Optimization using MATLAB
Solve this system for x, y, and z using the Gauss elimination method:
.
2.5 4.2 5.6 .. 26.50
3.4 4.7 2.8 ...
22.86
.
3.6 6.1 3.7 .. 29.12
..
2.5 4.2 5.6 . 26.50
.
≈ 0 −1.012 −4.816 .. −13.18
..
0 0.052 −4.364 . −9.04
Matrices and Linear Systems 193
.
2.5 4.2 5.6 .. 26.50
.
≈
0 −1.012 −4.816 .. −13.18
.
..
0 0 −4.612 . −9.717
Now use backward substitution to get the solution of the system:
For example, for the reaction in which hydrogen gas (H2 ) and oxygen
(O2 ) combine to form water (H2 O), a balanced chemical equation is
2H2 + O2 −→ 2H2 O,
Nitrogen: w = 2y
Hydrogen: 3w = 2z
Oxygen: 2x = z.
w − 2y = 0
3w − 2z = 0
2x − z = 0.
.
1 0 −2 0 .. 0
.
3 0
0 −2 .. 0
.
.
0 2 0 −1 .. 0
Solve this system for w, x, y, and z using the Gauss elimination method
with partial pivoting:
.. .
3 0 0 −2 . 0 3 0 0 −2 .. 0
2 .. .
0 0 −2
3 ≈ 0 2
. 0 0 −1 .. 0 .
. 2 ..
0 2 0 −1 .. 0 0 0 −2 3
. 0
Matrices and Linear Systems 195
Now use backward substitution to get the solution of the homogeneous sys-
tem:
−2y + 32 z = 0 gives y = 31 z
2x − z=0 gives x = 12 z
3w − 2z = 0 gives w = 32 z.
The smallest positive value of z that will produce integer values for all
four variables is the least common denominator of the fractions 23 , 12 , and
1
3
—namely, 6—which gives
w = 4, x = 3, y = 2, z = 6.
Therefore,
4N H3 + 3O2 −→ 2N2 + 6H2 O
is the balanced chemical equation. •
Solve this system for x, y, and z using the Gauss elimination method:
..
15 12 10 . 1250
..
15 12 10 . 1250
..
0 13 1 .
. 200
4 4.5 3 . 400 ≈
.
10 3 3
..
5 2.5 2.5 . 320
3 5 .. 290
0 − − . −
2 6 3
.
15 12 10 .. 1250
13 1 .
.. 200
≈ 0 .
10 3 3
35 .. 770
0 0 − . −
78 39
Now use backward substitution to get the solution of the system:
35 770
− z=− gives z = 44
78 39
13 1 200
y+ z= gives y = 40
10 3 3
Example 1.65 (Weather) The average of the temperature for the cities
of Jeddah, Makkah, and Riyadh was 50o C during a given summer day. The
temperature in Makkah was 5o C higher than the average of the temperatures
of the other two cities. The temperature in Riyadh was 5o C lower than the
average temperature of the other two cities. What was the temperature in
each of the cities?
50o C. On the other hand, the temperature in Makkah exceeds the average
temperature of Jeddah and Riyadh, (x+z)
2
, by 5o C. So, y = (x+z)
2+5
. Likewise,
(x+y)
we have z = 2−5 . So, the system becomes
(x + y + z)
= 50
3
(x + z)
y = +5
2
(x + y)
z = − 5.
2
Rewriting the above system in standard form, we get
x + y + z = 150
−x + 2y − z = 10
−x − y + 2z = −10.
Solve this system for x, y, and z using the Gauss elimination method:
.. ..
1 1 1 . 150 1 1 1 . 150
.. .
−1
2 −1 . 10 ≈ 0 3
0 .. 160 .
.. .
−1 −1 2 . −10 0 0 3 .. 140
Thus, the temperature in Jeddah was 50o C and the temperatures in Makkah
and Riyadh were approximately, 53o C and 470 C, respectively. •
following rates: the dollar was 60 rupees, 0.6 pounds, and 3.75 riyals. The
second time he exchanged a total of $25500 at these rates: the dollar was
65 rupees, 0.56 pounds, and 3.76 riyals. The third time he exchanged again
a total of $25500 at these rates: the dollar was 65 rupees, 0.6 pounds, and
3.75 riyals. How many rupees, pounds, and riyals did he buy each time?
1 1 1
x+ y+ z = 26000.
60 0.8 3.75
The same reasoning applies to the other two purchases, and we get the
system
1 5 4
x + y + z = 26000
60 3 15
1 25 25
x + y + z = 25500
65 14 94
1 5 4
x + y + z = 25500.
65 3 15
Solve this system for x, y, and z using the Gauss elimination method:
1 5 4 .. 1 5 4 ..
60 . 26000 . 26000
3 15 60 3 15
1 25 25 .. 45 121 ..
65 . 25500 ≈ 0 . 1500
14 94
182 6110
1 5 4 .. 5 4 ..
. 25500 0 . 1500
65 3 15 39 195
Matrices and Linear Systems 199
1 5 4 ..
60 . 26000
3 15
45 121 ..
≈ 0
. 1500 .
182 6110
13 .. 6500
0 0 .
1269 9
Now use backward substitution to get the solution of the system:
13 6500
z= gives z = 70500
1269 9
45 121
y+ z = 1500 gives y = 420
182 6110
1 5 4
x + y + z = 26000 gives x = 390000.
60 3 15
Therefore, each time he bought 390000 rupees, 420 pounds, and 70500 riyals
for his trips. •
Example 1.67 (Inheritance) A father plans to distribute his estate, worth
SR234,000, between his four daughters as follows: 23 of the estate is to be
split equally among the daughters. For the rest, each daughter is to receive
SR3,000 for each year that remains until her 21st birthday. Given that
the daughters are all 3 years apart, how much would each receive from her
father’s estate? How old are the daughters now?
Solve this system for x1 , x2 , x3 , and x4 using the Gauss elimination method
with partial pivoting:
.. ..
1 1 1 1 . 78, 000 1
..
1 1 1 . 78, 000
..
0 0 −1 1 . 9, 000 0 0 −1 1 . 9, 000
. ≈ ..
1 0 .. 9, 000
0 −1
0 −1 1 0 . 9, 000
.. ..
−1 1 0 0 . 9, 000 0 2 1 1 . 87, 000
and
..
..
1 1 1 1 . 78, 000 1 1 1 1 . 78, 000
.
.. 0 2 1 1 .. 87, 000
0 2 1 1 . 87, 000
3 1 .. ≈ 0
3 1 ..
.
0 0 . 52, 500 0 . 52, 500
2 2
2 2
.. 4 ..
0 0 −1 1 . 9, 000 0 0 0 . 44, 000
3
Now use backward substitution to get the solution of the system:
4
w = 44, 000 gives w = 33, 000
3
3 1
z + w = 52, 500 gives z = 24, 000
2 2
18x + 9y + 9z + 12w = 72
6x + 6y + 9z + 12w = 45
6x + 12y + 6z + 9w = 42
6x + 9y + 9z + 18w = 60.
Solve this system for x, y, z, and w using the Gauss elimination method:
.. ..
18 9 9 12 .
..
72 18 9 9 12 .
..
72
6 6 9 12 . 45 0 3 6 8 . 21
.. ≈ ..
6 12 6 9 . 42 0 9 3 5 . 18
.. ..
6 9 9 18 . 60 0 6 6 14 . 36
202 Applied Linear Algebra and Optimization using MATLAB
and
..
..
18 9 9 12 . 72
18 9 9 12
..
. 72
..
0 3
6 8 . 21
0 3 6 8 . 21
≈ .
.. ..
0 0 −15 −19 . −45 0 0 −15 −19 . −45
28 ..
..
0 0 −6 −2 . 36 0 0 0 . 12
5
Now use backward substitution to get the solution of the system:
28 15
w = 12 gives w =
5 7
2
−15z − 19w = −45 gives z =
7
5
3y + 6z + 8w = 21 gives y =
7
29
18x + 9y + 9z + 12w = 72 gives x = .
14
29
Thus, the amount in ounces of foods A, B, C, and D are x = 14
,y = 75 , z =
2
7
, and w = 15
7
, respectively. •
1.7 Summary
The basic methods for solving systems of linear algebraic equations were
discussed in this chapter. Since these methods use matrices and determi-
nants, the basic properties of matrices and determinants were presented.
Several direct solution methods were also discussed. Among them were
Cramer’s rule, Gaussian elimination and its variants, the Gauss–Jordan
method, and the LU decomposition method. Cramer’s rule is impracti-
cal for solving systems with more than three or four equations. Gaussian
elimination is the best choice for solving linear systems. For systems of
equations having a constant coefficients matrix but many right-hand side
Matrices and Linear Systems 203
1.8 Problems
1. Determine the matrix C given by the following expression
C = 2A − 3B,
4. Let
1 2 3 1 1 2 1 0 1
A = 0 −1 2 , B = −1 1 −1 , C = 0 1 2 .
2 0 2 1 0 2 2 0 1
6. Find the values of a and b such that each of the following matrices is
symmetric:
1 3 5 −2 a + b 2
(a) A = a + 2 5 6 , (b) B = 3 4 2a + b ,
b+1 6 7 2 5 −3
1 4 a−b 1 a − 4b 2
(c) C = 4 2 a + 3b , (d) D = 2 8 6 .
7 3 4 7 a − 7b 8
(a)
1 −5 0 −4
A= , B= ,
5 0 4 0
(b)
1 9 1 6
C= , D= ,
−9 7 −6 2
(c)
0 2 −2 3 −3 −3
E = −2 0 4 , F = 3 3 −3 ,
2 −4 0 3 3 3
(d)
1 −5 1 2 8 6
G= 5 1 4 , H = −8 4 2 .
−1 −4 1 −6 −2 5
(a)
1 0 0 1 0 8
A = 0 1 0 , B = 0 1 2 ,
0 0 3 0 0 0
206 Applied Linear Algebra and Optimization using MATLAB
(b)
1 2 3 0 1 2 0 0 1
C = 0 0 0 1 , D = 0 0 1 0 1 ,
0 0 0 1 0 0 0 1 0
(c)
1 4 5 6
0 1 0 0 0 3
1 7 8
E=
0
, F = 0 0 1 0 4 ,
0 1 9
0 0 0 1 5
0 0 0 0
(d)
1 0 0 3 0 0 0 0 0
0 1 0 4 0 0 1 2 4
G=
0
, H= .
0 0 5 0 0 0 1 0
0 0 0 6 0 0 0 0 0
9. Find the row echelon form of each of the following matrices using
elementary row operations, and then solve the linear system:
(a)
0 1 2 1
A = 2 3 4 , b = −1 .
1 3 2 2
(b)
1 2 3 1
A = 0 3 1 , b = 0 .
−1 4 5 −3
(c)
0 −1 0 1
A= 3 0 1 , b = 3 .
0 1 1 2
(d)
0 −1 2 4 2
2 3 5 6 1
A= , b=
−1 .
1 3 −2 4
1 2 −1 3 2
Matrices and Linear Systems 207
10. Find the row echelon form of each of the following matrices using
elementary row operations, and then solve the linear system:
(a)
1 4 −2 5
A = 2 3 2 , b = −3 .
6 4 1 4
(b)
2 2 7 3
A = 0 3 2 , b = 2 .
3 2 1 5
(c)
0 −1 0 1
A= 5 0 2 , b = 1 .
−1 1 4 1
(d)
1 1 2 4 11
1 3 4 5 7
A=
1
, b=
6 .
4 2 4
2 2 −1 3 4
11. Find the reduced row echelon form of each of the following matrices
using elementary row operations, and then solve the linear system:
(a)
1 2 3 4
A = −1 2 1 , b = 3 .
0 1 2 1
(b)
0 1 4 1
A = 2 1 −1 , b = 1 .
1 3 4 −1
208 Applied Linear Algebra and Optimization using MATLAB
(c)
0 −1 3 2 6
3 2 5 4 4
A=
−1
, b=
4 .
3 1 2
2 3 4 1 4
(d)
1 2 −4 1 1
−2 0 2 3 −1
A= , b=
2 .
0 1 −1 2
2 3 0 −1 4
14. Let
1 1 1 0
A= , B= ,
0 1 1 1
then show that (AB)−1 = B −1 A−1 .
15. Evaluate the determinant of each of the following matrices using the
Gauss elimination method:
3 1 −1 4 1 6 17 46 7
A= 2 0 4 , B = −3 6 4 , C = 20 49 8 .
1 −5 1 5 0 9 23 52 19
Matrices and Linear Systems 209
16. Evaluate the determinant of each of the following matrices using the
Gauss elimination method:
4 2 5 −1 4 −2 5 −3
2 5 4 6
, B = 1 8 12 7
A= 4 5 1
,
3 1 4 3 6
11 7 1 1 5 3 −3 6
13 22 −12 8 9 11 2 8
15 10 33 4 15 1 3 12
C=
9 −12
, D=
9 −12 5 17 .
5 7
15 33 −19 26 13 17 21 15
17. Find all zeros (values of x such that f (x) = 0) of polynomial f (x) =
det(A), where
x−1 3 2
A= 3 x 1 .
2 1 x−2
18. Find all zeros (values of x such that f (x) = 0) of polynomial f (x) =
det(A), where
x 0 1
A = 2 1 3 .
0 x 2
19. Find all zeros (values of x such that f (x) = 0) of polynomial f (x) =
det(A), where
x −8 5 2
−3 x 2 1
A= 3
.
4 x 1
3 6 −5 17
23. Find the inverse and determinant of the adjoint matrix of each of the
following matrices:
4 1 5 3 4 −2 1 2 4
A = 5 6 3 , B = 2 5 4 , C = 1 4 0 .
5 4 4 7 −3 4 3 1 1
24. Find the inverse and determinant of the adjoint matrix of each of the
following matrices:
3 2 5 5 3 −2 1 2 3
A = 2 5 4 , B = 3 5 6 , C = 4 5 6 .
5 4 6 −2 6 5 7 8 8
25. Find the inverse of each of the following matrices using the determi-
nant:
0 4 2 −4
0 1 5 2 4 −2 6 1 4 −3
A = 3 1 2 , B = −4 7 5 , C = 4
.
3 1 3
2 3 4 5 −4 4
8 4 −3 2
(a)
x1 − 2x2 + x3 = 0
x1 + x2 + 3x3 = 0
2x1 + 3x2 − 5x3 = 0.
212 Applied Linear Algebra and Optimization using MATLAB
(b)
x1 − 5x2 + 3x3 = 0
2x1 + 3x2 + 2x3 = 0
x1 − 2x2 − 4x3 = 0.
(c)
3x1 + 4x2 − 2x3 = 0
2x1 − 5x2 − 4x3 = 0
3x1 − 2x2 + 3x3 = 0.
(d)
x1 + x2 + 3x3 − 2x4 = 0
x1 + 2x2 + 5x3 + x4 = 0
x1 − 3x2 + x3 + 2x4 = 0.
27. Find value(s) of α such that each of the following homogeneous linear
systems has a nontrivial solution:
(a)
2x1 − (1 − 3α)x2 = 0
x1 + αx2 = .
(b)
2x1 + 2αx2 − x3 = 0
x1 − 2x2 + x3 = 0
αx1 + 2x2 − 3x3 = 0.
(c)
x1 + 2x2 + 4x3 = 0
3x1 + 7x2 + αx3 = 0
3x1 + 3x2 + 15x3 = 0.
(d)
x1 + x2 + 2x3 − 3x4 = 0
x1 + 2x2 + x3 − 2x4 = 0
3x1 + x2 + αx3 + 3x4 = 0
3x1 + x2 + αx3 + 3x4 = 0
2x1 + 3x2 + x3 + αx4 = 0.
Matrices and Linear Systems 213
28. Using the matrices in Problem 15, solve the following systems using
the matrix inversion method:
29. Solve the following systems using the matrix inversion method:
(a)
x1 + 3x2 − x3 = 4
5x1 − 2x2 − x3 = −2
2x1 + 2x2 + x3 = 9.
(b)
x1 + x2 + 3x3 = 2
5x1 + 3x2 + x3 = 3
2x1 + 3x2 + x3 = −1.
(c)
4x1 + x2 − 3x3 = −1
3x1 + 2x2 − 6x3 = −2
x1 − 5x2 + 3x3 = −3.
(d)
7x1 + 11x2 − 15x3 = 21
3x1 + 22x2 − 18x3 = 12
2x1 − 13x2 + 9x3 = 16.
30. Solve the following systems using the matrix inversion method:
(a)
3x1 − 2x2 − 4x3 = 7
5x1 − 2x2 − 3x3 = 8
7x1 + 4x2 + 2x3 = 9.
(b)
−3x1 + 4x2 + 3x3 = 11
5x1 + 3x2 + x3 = 12
x1 + x2 + 5x3 = 10.
214 Applied Linear Algebra and Optimization using MATLAB
(c)
x1 + 42 − 8x3 = 7
2x1 + 7x2 − 5x3 = −5
3x1 − 6x2 + 6x3 = 4.
(d)
17x1 + 18x2 − 19x3 = 10
43x1 + 22x2 − 14x3 = 11
25x1 − 33x2 + 21x3 = 12.
31. Solve the following systems using the matrix inversion method:
(a)
2x1 + 3x2 − 4x3 + 4x4 = 11
x1 + 3x2 − 4x3 + 2x4 = 12
4x1 + 3x2 + 2x3 + 3x4 = 14
3x1 − 4x2 + 5x3 + 6x4 = 15.
(b)
7x1 + 13x2 + 12x3 + 9x4 = 21
3x1 + 23x2 − 5x3 + 2x4 = 10
4x1 − 7x2 + 22x3 + 3x4 = 11
3x1 − 4x2 + 25x3 + 16x4 = 10.
(c)
12x1 + 6x2 + 5x3 − 2x4 = 21
11x1 + 13x2 + 7x3 + 2x4 = 22
14x1 + 9x2 + 2x3 − 6x4 = 23
7x1 − 24x2 − 7x3 + 8x4 = 24.
(d)
15x1 − 26x2 + 15x3 − 11x4 = 17
14x1 + 15x2 + 7x3 + 7x4 = 18
17x1 + 14x2 − 22x3 − 16x4 = 19
21x1 − 12x2 − 7x3 + 8x4 = 20.
(a)
3x1 + 4x2 + 5x3 = 1
3x1 + 2x2 + x3 = 2
4x1 + 3x2 + 5x3 = 3.
(b)
x1 − 4x2 + 2x3 = 4
−4x1 + 5x2 + 6x3 = 0
7x1 − 3x2 + 5x3 = 4.
(c)
6x1 + 7x2 + 8x3 = 1
−5x1 + 3x2 + 2x3 = 1
x1 + 2x2 + 3x3 = 1.
(d)
x1 + 3x2 − 4x3 + 5x4 = 2
6x1 − x2 + 6x3 + 3x4 = −3
2x1 + x2 + 3x3 + 2x4 = 4
x1 + 5x2 + 6x3 + 7x4 = 2.
(a)
2x1 − 2x2 + 8x3 = 1
5x1 + 6x2 + 5x3 = 2
7x1 + 7x2 + 9x3 = 3.
(b)
3x1 − 3x2 + 12x3 = 14
−4x1 + 5x2 + 16x3 = 18
x1 − 15x2 + 24x3 = 19.
216 Applied Linear Algebra and Optimization using MATLAB
(c)
9x1 − 11x2 + 12x3 = 3
−5x1 + 3x2 + 2x3 = 4
7x1 − 12x2 + 13x3 = 5.
(d)
11x1 + 3x2 − 13x3 + 15x4 = 22
26x1 − 5x2 + 6x3 + 13x4 = 23
22x1 + 6x2 + 13x3 + 12x4 = 24
17x1 − 25x2 + 16x3 + 27x4 = 25.
36. Use the simple Gaussian elimination method to show that the fol-
lowing system does not have a solution:
3x1 + x2 = 1.5
2x1 − x2 − x3 = 2
4x1 + 3x2 + x3 = 0.
38. Solve the following systems using the simple Gaussian elimination
method:
(a)
x1 − x2 = −2
−x1 + 2x2 − x3 = 5
4x1 − x2 + 4x3 = 1.
(b)
3x1 + x2 − x3 = 5
5x1 − 3x2 + 2x3 = 7
2x1 − x2 + x3 = 3.
(c)
3x1 + x2 + x3 = 2
2x1 + 2x2 + 4x3 = 3
4x1 + 9x2 + 16x3 = 1.
Matrices and Linear Systems 217
(d)
2x1 + x2 + x3 − x4 = 9
x1 + 9x2 + 8x3 + 4x4 = 11
−x1 + 3x2 + 5x3 + 2x4 = 10
5x1 + x2 + x4 = 12.
39. Solve the following systems using the simple Gaussian elimination
method:
(a)
2x1 + 5x2 − 4x3 = 3
2x1 + 2x2 − x3 = 1
3x1 + 2x2 − 3x3 = −5.
(b)
2x2 − x3 = 1
3x1 − x2 + 2x3 = 4
x1 + 3x2 − 5x3 = 1.
(c)
x1 + 2x2 = 3
−x1 − 2x3 = −5
−3x1 − 5x2 + x3 = −4.
(d)
3x1 + 2x2 + 4x3 − x4 = 2
x1 + 4x2 + 5x3 + x4 = 1
4x1 + 5x2 + 4x3 + 3x4 = 5
2x1 + 3x2 + 2x3 + 4x4 = 6.
40. For what values of a and b does the following linear system have no
solution or infinitely many solutions:
(a)
2x1 + x2 + x3 = 2
−2x1 + x2 + 3x3 = a
2x1 − x3 = b.
218 Applied Linear Algebra and Optimization using MATLAB
(b)
2x1 + 3x2 − x3 = 1
x1 − x2 + 3x3 = a
3x1 + 7x2 − 5x3 = b.
(c)
2x1 − x2 + 3x3 = 3
3x1 + x2 − 5x3 = a
−5x1 − 5x2 + 21x3 = b.
(d)
2x1 − x2 + 3x3 = 5
4x1 + 2x2 + bx3 = 6
−2x1 + ax2 + 3x3 = 4.
41. Find the value(s) of α so that each of the following linear systems
has a nontrivial solution:
(a)
2x1 + 2x2 + 3x3 = 1
3x1 + αx2 + 5x3 = 3
x1 + 7x2 + 3x3 = 2.
(b)
x1 + 2x2 + x3 = 2
x1 + 3x2 + 6x3 = 5
2x1 + 3x2 + αx3 = 6.
(c)
αx1 + x2 + x3 = 7
x1 + x2 − x3 = 2
x1 + x2 + αx3 = 1.
(d)
2x1 + αx2 + 3x3 = 9
3x1 − 4x2 − 5x3 = 11
4x1 + 5x2 + αx3 = 12.
Matrices and Linear Systems 219
42. Find the inverse of each of the following matrices by using the simple
Gauss elimination method:
3 3 3 5 3 2 1 2 3
A = 0 2 2 , B = 3 2 2 , C = 2 5 2 .
2 4 5 2 6 5 3 4 3
43. Find the inverse of each of the following matrices by using the simple
Gauss elimination method:
3 2 3 1 −3 2 5 2 3
A = 4 2 2 , B = 3 2 6 , C = 2 5 5 .
2 4 3 2 −6 5 3 2 4
48. Solve the following linear systems using Gaussian elimination with
partial and without pivoting:
(a)
1.001x1 + 1.5x2 = 0
2x1 + 3x2 = 1.
220 Applied Linear Algebra and Optimization using MATLAB
(b)
x1 + 1.001x2 = 2.001
x1 + x2 = 2.
(c)
6.122x1 + 1500.5x2 = 1506.622
2000x1 + 3x2 = 2003.
49. The elements of matrix A, the Hilbert matrix, are defined by
aij = 1/(i + j − 1), for i, j = 1, 2, . . . , n.
Find the solution of the system Ax = b for n = 4 and b = [1, 2, 3, 4]T
using Gaussian elimination by partial pivoting.
50. Solve the following systems using the Gauss–Jordan method:
(a)
x1 + 4x2 + x3 = 1
2x1 + 4x2 + x3 = 9
3x1 + 5x2 − 2x3 = 11.
(b)
x1 + x 2 + x3 = 1
2x1 − x2 + 3x3 = 4
3x1 + 2x2 − 2x3 = −2.
(c)
2x1 + 3x2 + 6x3 + x4 = 2
x1 + x2 − 2x3 + 4x4 = 1
3x1 + 5x2 − 2x3 + 2x4 = 11
2x1 + 2x2 + 2x3 − 3x4 = 2.
51. The following sets of linear equations have a common coefficients ma-
trix but different right-side terms:
(a)
2x1 + 3x2 + 5x3 = 0
3x1 + x2 − 2x3 = −2
x1 + 3x2 + 4x3 = −3.
Matrices and Linear Systems 221
(b)
2x1 + 3x2 + 5x3 = 1
3x1 + x2 − 2x3 = 2
x1 + 3x2 + 4x3 = 4.
(c)
2x1 + 3x2 + 5x3 = −5
3x1 + x2 − 2x3 = 6
x1 + 3x2 + 4x3 = −1.
The coefficients and the three sets of right-side terms may be com-
bined into an augmented matrix of the form
..
2 3 5 . 0 1 −5
3 1 −2 ..
.
. −2 2 6
..
1 3 4 . −3 4 −1
If we apply the Gauss–Jordan method to this augmented matrix form
and reduce the first three columns to the unity matrix form, the solu-
tion for the three problems are automatically obtained in the fourth,
fifth, and sixth columns when elimination is completed. Calculate
the solution in this way.
52. Calculate the inverse of each matrix using the Gauss–Jordan method:
5 −2 0 0
3 −9 5 1 4 5 −2 5 −2 0
(a) 0 5 1 , (b) 2 1 2 , (c) 0 −2
.
5 −2
−1 6 3 8 1 1
0 0 −2 5
53. Find the inverse of the Hilbert matrix of size 4 × 4 using the Gauss–
Jordan method. Then solve the linear system Ax = [1, 2, 3, 4]T .
54. Find the LU decomposition of each matrix A using Doolittle’s method
and then solve the systems:
(a)
2 −1 1 4
A= −3 4 −1 , b= 5 .
1 −1 1 6
222 Applied Linear Algebra and Optimization using MATLAB
(b)
7 6 5 2
A = 5 4 3 , b = 1 .
3 7 6 2
(c)
2 2 2 0
A = 1 2 1 , b = −4 .
3 3 4 1
(d)
2 4 −6 −4
A= 1 5 3 , b = 10 .
1 3 2 5
(e)
1 −1 0 2
A = 2 −1 1 , b = 4 .
2 −2 −1 3
(f )
1 5 3 4
A = 2 4 6 , b = 11 .
1 3 2 5
55. Find the LU decomposition of each matrix A using Doolittle’s method,
and then solve the systems:
(a)
3 −2 1 1 3
−3 7 4 −3 2
A=
2 −5 3
, b=
1 .
4
7 −3 2 4 2
(b)
2 −4 5 3 6
3 5 −4 3 5
A=
1
, b=
2 .
6 2 6
7 2 5 1 4
Matrices and Linear Systems 223
(c)
2 2 3 −2 10
10 2 13 11 14
A=
2
, b=
11 .
5 4 6
1 −4 −2 7 9
(d)
5 12 4 −11 44
21 15 13 23 33
A=
31
, b=
55 .
33 12 22
−17 15 14 11 22
(e)
1 −1 10 8 −2
12 −17 11 22 7
A= , b=
6 .
22 31 13 −1
8 24 13 9 5
(f )
41
41 25 23 −18 1
A = 2 13 −16 12 , 15 .
b=
11 13 9 7
13
1 −1 2
(a) A = −1 3 −1 .
α −2 3
1 5 7
(b) A = 4 4 α .
−2 α 9
2 −4 α
(c) A= 2 4 3 .
4 −2 5
2 α 1−α
(d) A= 2 5 −2 .
2 5 4
1 −1 3
(e) A= 3 2 3 .
4 α−2 7
1 5 α
(f ) A= 1 4 α − 2 .
1 −2 8
58. Use the smallest positive integer to find the unique solution of each of
the linear systems of Problem 56 using LU decomposition by Doolit-
tle’s method:
(a)
3 4 3 4 −2 3
A = 2 3 3 , B= 5 2 −3 .
1 3 5 4 3 6
(b)
2 5 4 3 2 −6
A = 2 1 6 , B = 2 2 −5 .
3 2 7 3 4 7
(c)
1 −5 4 4 7 −6
A= 2 3 −4 , B= 5 5 −5 .
3 2 6 6 −4 9
(d)
2 3 4 5
3 −1 4 3 1 2 4
A= 2 2 −1 , B=
3
.
1 1 1
3 2 2
4 3 1 2
226 Applied Linear Algebra and Optimization using MATLAB
(a)
2 3 4
A = 3 5 2 .
4 2 6
(b)
3 −2 4
A = −2 2 1 .
4 1 3
(c)
2 1 −1
A= 1 3 2 .
−1 2 2
(d)
1 −2 3 4
−2 3 4 5
A= .
3 4 5 −6
4 5 −6 7
(a)
1 −1 1 2
A = −1 5 −1 , b = 2 .
1 −1 10 2
(b)
10 2 1 7
A = 2 10 3 , b = −4 .
1 3 10 3
(c)
4 2 3 1
A = 2 17 1 , b = 2 .
3 1 5 5
(d)
3 4 −6 0 4
4 5 3 1 5
A=
−6
, b=
2 .
3 3 1
0 1 1 3 3
(a)
5 −1 1 5
A = −3 5 −2 , b = 7 .
2 −1 7 9
(b)
6 2 −3 5
A = 3 12 −4 , b = 2 .
6 3 13 4
228 Applied Linear Algebra and Optimization using MATLAB
(c)
5 2 −5 3
A= 2 4 4 , b = 11 .
−3 −2 7 14
(d)
1 4 −6 0 12
2 2 3 3 13
A=
−3
, b=
14 .
6 7 1
0 2 −3 5 15
(a)
3 −1 0 1
A = −1 3 −1 , b = 1 .
0 −1 3 1
(b)
2 3 0 0 6
3 2 3 0 7
A=
0
, b=
5 .
3 2 3
0 0 3 2 3
(c)
4 −1 0 0 1
−1 4 −1 0 1
A=
0 −1
, b=
1 .
4 −1
0 0 −1 4 1
(d)
2 3 0 0 1
3 5 4 0 2
A=
0
, b=
3 .
4 6 3
0 0 3 4 4
(a)
4 −2 0 5
A = −2 5 −2 , b = 6 .
0 −2 6 7
(b)
8 1 0 0 2
1 8 1 0 2
A=
0
, b=
2 .
1 8 1
0 0 1 8 2
(c)
5 −3 0 0 7
−3 6 −2 0 −5
A=
0 −2
, b=
4 .
7 −5
0 0 −5 8 2
(d)
2 −4 0 0 11
−4 5 7 0 12
A=
0
, b=
13 .
7 6 2
0 0 2 8 14
67. Find kxk1 , kxk2 , and kxk∞ for the following vectors:
(a)
[2, −1, −6, 3]T .
(b)
[sin k, cos k, 3k ]T , for a fixed integer k.
68. Find k.k1 , k.k∞ and [Link] for the following matrices:
3 1 −1 4 1 6
A= 2 0 4 , B = −3 6 4 ,
1 −5 1 5 0 9
3 11 −5 2
17 46 7 6 8 −11 6
C= 20 49 8 , D = −4 −8
.
10 14
23 52 9
13 14 −12 9
230 Applied Linear Algebra and Optimization using MATLAB
(a)
0.89x1 + 0.53x2 = 0.36
0.47x1 + 0.28x2 = 0.19
Matrices and Linear Systems 231
x = [1, −1]T
x∗ = [0.702, −0.500]T
(b)
0.986x1 + 0.579x2 = 0.235
0.409x1 + 0.237x2 = 0.107
x = [2, −3]T
x∗ = [2.110, −3.170]T
(c)
1.003x1 + 58.090x2 = 68.12
5.550x1 + 321.8x2 = 377.3
x = [10, 1]T
x∗ = [−10, 1]T
1.01x1 + 0.99x2 = 2
0.99x1 + 1.01x2 = 2.
1 kA − Bk
≤ .
K(A) kAk
(a)
K(A) ≥ 1 and K(B) ≥ 1.
(b)
K(AB) ≤ K(A)K(B).
80. If kAk < 1, then show that the matrix (I − A) is nonsingular and
1
k(I − A)−1 k ≤ .
1 − kAk
85. Determine the currents through the various branches of the electrical
network in Figure 1.8:
Note how the current through the branch AB is reversed in (b). What
would the voltage of C have to be for no current to pass through AB?
Matrices and Linear Systems 235
87. Figure 1.10 represents the traffic entering and leaving a “roundabout”
road junction. Such junctions are very common in Europe. Construct
a mathematical model that describes the flow of traffic along the vari-
ous branches. What is the minimum flow theoretically possible along
the branch BC? Is this flow ever likely to be realized in practice?
(a) C7 H6 O2 + O2 −→ H2 O + CO2 .
94. The average of the temperature for the cities of Jeddah, Makkah,
and Riyadh was 15o C during a given winter day. The temperature
in Makkah was 6o C higher than the average of the temperatures
of the other two cities. The temperature in Riyadh was 6o C lower
than the average temperature of the other two cities. What was the
temperature in each one of the cities?
97. A biologist has placed three strains of bacteria (denoted by I, II, and
III) in a test tube, where they will feed on three different food sources
(A, B, and C). Each day 2300 units of A, 800 units of B, and 1500
units of C are placed in the test tube, and each bacterium consumes
a certain number of units of each food per day, as shown in the given
table. How many bacteria of each strain can coexist in the test tube
and consume all the food?
98. Al-karim hires three types of laborers, I, II, and III, and pays them
SR20, SR15, and SR10 per hour, respectively. If the total amount
paid is SR20,000 for a total of 300 hours of work, find the possible
number of hours put in by the three categories of workers if the
category III workers must put in the maximum amount of hours.
Chapter 2
2.1 Introduction
The methods discussed in Chapter 1 for the solution of the system of linear
equations have been direct, which required a finite number of arithmetic
operations. The elimination methods for solving such systems usually
yield sufficiently accurate solutions for approximately 20 to 25 simulta-
neous equations, where most of the unknowns are present in all of the
equations. When the coefficients matrix is sparse (has many zeros), a con-
siderably large number of equations can be handled by the elimination
methods. But these methods are generally impractical when many hun-
dreds or thousands of equations must be solved simultaneously.
There are, however, several methods that can be used to solve large
numbers of simultaneous equations. These methods, called iterative meth-
ods, are methods by which an approximation to the solution of a system
243
244 Applied Linear Algebra and Optimization using MATLAB
of linear equations may be obtained. The iterative methods are used most
often for large, sparse systems of linear equations and they are efficient in
terms of computer storage and time requirements. Systems of this type
arise frequently in the numerical solutions of boundary value problems
and partial differential equations. Unlike the direct methods, the iterative
methods may not always yield a solution, even if the determinant of the
coefficients matrix is not zero.
Ax = b (2.1)
x = Tx + c (2.2)
for some square matrix T and vector c. After the initial vector x(0) is se-
lected, the sequence of approximate solution vectors is generated by com-
puting
x(k+1) = T x(k) + c, for k = 0, 1, 2, . . . . (2.3)
Among them, the most useful methods are the Jacobi method, the
Gauss–Seidel method, the Successive Over-Relaxation (SOR) method, and
the conjugate gradient method.
A = L + D + U, (2.5)
Iterative Methods for Linear Systems 245
or in matrix form
Dx = b − (L + U )x.
Divide both sides of the above three equations by their diagonal elements,
a11 , a22 , and a33 , respectively, to get
1 h i
x1 = b1 − a12 x2 − a13 x3
a11
1 h i
x2 = b2 − a21 x1 − a23 x3
a22
1 h i
x3 = b3 − a31 x1 − a32 x2 ,
a33
x = D−1 [b − (L + U )x].
h iT
(0) (0) (0) (0)
Let x = x1 , x2 , x3
be an initial solution of the exact solution x of
the linear system (2.1). Then define an iterative sequence
or in matrix form
where k is the number of iterative steps. Then the form (2.7) is called
the Jacobi formula for the system of three equations and (2.8) is called its
matrix form. For a general system of n linear equations, the Jacobi method
Iterative Methods for Linear Systems 247
is defined by
" i−1 n
#
(k+1) 1 X (k)
X (k)
xi = bi − aij xj − aij xj (2.9)
aii j=1 j=i+1
i = 1, 2, . . . , n, k = 0, 1, 2, . . .
or
x1 c1 0 −t12 · · · −t1n x1
x2 −t21
c2 0 · · · −t2n x2
= + .. .. ,
.. .. .. .. ..
. . . . . . .
xn (k+1) cn −tn1 −tn2 · · · 0 xn (k)
(2.11)
where the Jacobi iteration matrix TJ and vector c are defined as follows:
tij = 0, i=j
bi
ci = , i = 1, 2, . . . , n.
aii
The Jacobi iterative method is sometimes called the method of simultane-
ous iterations, because all values of xi are iterated simultaneously. That
(k+1) (k)
is, all values of xi depend only on the values of xi .
248 Applied Linear Algebra and Optimization using MATLAB
Note that the diagonal elements of the Jacobi iteration matrix TJ are
(0)
always zero. As usual with iterative methods, an initial approximation xi
must be supplied. If we don’t have knowledge of the exact solution, it is
(0)
conventional to start with xi = 0, for all i. The iterations defined by
(2.9) are stopped when
kx(k+1) − x(k) k
< (2.14)
kx(k+1) k
Example 2.1 Solve the following system of equations using the Jacobi it-
erative method using = 10−5 in the l∞ -norm:
Note that the Jacobi method converged and after 15 iterations we ob-
tained the good approximation [2.24373, 2.93123, 3.57829, 4.18940]T of
the given system having the exact solution [2.24374, 2.93124, 3.57830,
4.18941]T . Ideally, the iterations should stop automatically when we ob-
tain the required accuracy using one of the stopping criteria mentioned in
(2.13) or (2.14). •
Example 2.2 Solve the following system of equations using the Jacobi it-
erative method:
Solution. Results for this linear system are listed in Table 2.2. Note that
in this case the Jacobi method diverges rapidly. Although the given linear
system is the same as the linear system of Example 2.1, the first and second
equations are interchanged. From this example we conclude that the Jacobi
iterative method is not always convergent.
Iterative Methods for Linear Systems 251
Program 2.1
MATLAB m-file for the Jacobi Iterative Method
function x=JacobiM(Ab,x,acc) % Ab = [A b]
[n,t]=size(Ab); b=Ab(1:n,t); R=1; k=1;
d(1,1:n+1)=[0 x]; while R > acc
for i=1:n
sum=0;
for j=1:n; if j ˜ =i
sum = sum + Ab(i, j) ∗ d(k, j + 1); end;
x(1, i) = (1/Ab(i, i)) ∗ (b(i, 1) − sum); end;end
k=k+1; d(k,1:n+1)=[k-1 x];
R=max(abs((d(k,2:n+1)-d(k-1,2:n+1))));
if k > 10 & R > 100
(‘Jacobi Method diverges’)
break; end; end; x=d;
bi
3. Compute the constant c = D−1 b = , for i = 1, 2, . . . , n.
aii
(k+1) (k)
5. Solve for the approximate solutions xi = TJ xi +c, i = 1, 2, . . . , n
k = 0, 1, . . ..
(k+1) (k)
6. Repeat step 5 until kxi − xi k < .
From the Jacobi iterative formula (2.9), it is seen that the new estimates
for solution x are computed from the old estimates and only when all
the new estimates have been determined are they then used in the right-
hand side of the equation to perform the next iteration. But the Gauss–
Seidel method is used to make use of the new estimates in the right-hand
side of the equation as soon as they become available. For example, the
Gauss–Seidel formula for the system of three equations can be defined as
an iterative sequence:
which are called the Gauss–Seidel iteration matrix and the vector, respec-
tively.
Example 2.3 Solve the following system of equations using the Gauss–
Seidel iterative method, with = 10−5 in the l∞ -norm:
Note that the Gauss–Seidel method converged for the given system and re-
quired nine iterations to obtain the approximate solution [2.24374, 2.93123,
3.57830, 4.18941]T , which is equal to the exact solution [2.24374, 2.93124,
3.57830, 4.18941]T up to six significant digits, which is six iterations less
than required by the Jacobi method for the same linear system. •
Example 2.4 Solve the following system of equations using the Gauss–
Seidel iterative method, with = 10−5 in the l∞ -norm:
Solution. Results for this linear system are listed in Table 2.4. Note that
in this case the Gauss–Seidel method diverges rapidly. Although the given
linear system is the same as the linear system of the previous Example 2.3,
the first and second equations are interchanged. From this example we
conclude that the Gauss–Seidel iterative method is not always convergent.•
256 Applied Linear Algebra and Optimization using MATLAB
Program 2.2
MATLAB m-file for the Gauss–Seidel Iterative Method
function x=GaussSM(Ab,x,acc) % Ab = [A b]
[n,t]=size(Ab); b=Ab(1:n,t);R=1; k=1;
d(1,1:n+1)=[0 x]; k=k+1; while R > acc
for i=1:n; sum=0; for j=1:n
if j <= i − 1; sum = sum + Ab(i, j) ∗ d(k, j + 1);
elseif j >= i + 1
sum = sum + Ab(i, j) ∗ d(k − 1, j + 1); end; end
x(1, i) = (1/Ab(i, i)) ∗ (b(i, 1) − sum);
d(k,1)=k-1; d(k,i+1)=x(1,i); end
R=max(abs((d(k,2:n+1)-d(k-1,2:n+1))));
k=k+1; if R > 100 & k > 10 (‘Gauss–Seidel method Diverges’)
break ;end;end;x=d;
Procedure 2.2 (Gauss–Seidel Method)
1. Check that the coefficient matrix A is strictly diagonally dominant
(for guaranteed convergence).
2. Initialize the first approximation x(0) ∈ R and preassigned accuracy
.
Iterative Methods for Linear Systems 257
The first and subsequent iterations are listed in Table 2.5. Now we solve
the same system by the Gauss–Seidel method and for the given system, the
Iterative Methods for Linear Systems 259
Gauss–Seidel formula is
The first and subsequent iterations are listed in Table 2.6. Note that the Ja-
cobi method diverged and the Gauss–Seidel method converged after 28 itera-
tions with the approximate solution [0.66150, −0.28350, 0.63177, 0.48758]T
of the given system, which has the exact solution [0.66169, −0.28358, 0.63184,
0.48756]T . •
Example 2.6 Solve the following system of equations using the Jacobi and
Gauss–Seidel iterative methods, using = 10−5 in the l∞ -norm and taking
the initial solution x(0) = [0, 0, 0, 0]T :
x1 + 2x2 − 2x3 = 1
x1 + x2 + x3 = 2
2x1 + 2x2 + x3 = 3
x1 + x2 + x3 + x 4 = 4.
Solution. First, we solve by the Jacobi method and for the given system,
the Jacobi formula is
The first and subsequent iterations are listed in Table 2.7. Now we solve
the same system by the Gauss–Seidel method and for the given system, the
262 Applied Linear Algebra and Optimization using MATLAB
Gauss–Seidel formula is
(k+1) 1h (k) (k)
i
x1 = 1 − 2x2 + 2x3
1
Jacobi method converged quickly (only five iterations) but the Gauss–Seidel
method diverged for the given system. •
0 0 0 0 2 0 6 0 0
A=L+U +D = 1 0 0 + 0 0 −2 + 0 7 0 .
3 −2 0 0 0 0 0 0 9
(a) Since the matrix form of the Jacobi iterative method can be written as
x(k+1) = TJ x(k) + c, k = 0, 1, 2, . . .
where
TJ = −D−1 (L + U ) and c = D−1 b
one can easily compute the Jacobi iteration matrix TJ and the vector c as
follows:
2 1
0 −6 0 6
1 2 2
TJ =
−7 0 and c =
7 .
7
3 2 1
− 0 −
9 9 9
Thus, the matrix form of the Jacobi iterative method is
2 1
0 − 0
6
6
1 2 (k)
2
x(k+1) =
−7 0 x + , k = 0, 1, 2.
7
7
3 2 1
− 0 −
9 9 9
(b) Now by writing the above iterative matrix form in component form, we
Iterative Methods for Linear Systems 265
have
1 1
0 − 0
x1
3 x1
6
x2 = − 1 2 2
0 x2 + ,
7 7
7
x3 1 2 x3 1
− 0 −
3 9 9
and it is equivalent to
1 1
x1 = − x2 +
3 6
1 2 2
x2 = − x1 + x3 +
7 7 7
1 2 1
x3 = − x1 + x2 − .
3 9 9
Now solving for x1 , x2 , and x3 , we get
1
x1 12
x2 = 1
,
4
x3 1
−
12
which is the exact solution of the given system.
(c) Since the error in the (n + 1)th step is defined as
e(k+1) = x − x(k+1) ,
we have
1 2 1
0 − 0
12
6
6
1 1
2 (k)
2
e(k+1) = − − 0 x + .
4 7
7
7
1 3 2 1
− − 0 −
12 9 9 9
266 Applied Linear Algebra and Optimization using MATLAB
1 2 1
0 − 0
12
6 12
1 1
2 1
e(k+1) = − − 0
4 7
7
4
1 3 2 1
− − 0 −
12 9 9 12
2
0 − 0
6
1 2
e(k)
+ − 0
7 7
3 2
− 0
9 9
1
6
2
+
7
1
−
9
or
2
0 − 0
6
1 2
e(k+1) =
−7 0 e(k) ,
7
3 2
− 0
9 9
(because x(k) = e(k) − x) which is the required error in the (n + 1)th step.
(d) Now finding the first approximation of the error, we have to compute
Iterative Methods for Linear Systems 267
the following:
2
0 − 0
6
1 2
e(1) =
−7 0 e(0) ,
7
3 2
− 0
9 9
where
e(0) = x − x(0) .
Using x(0) = [0, 0, 0]T , we have
1 1
12 0 12
(0)
1
− 0 = 1
e = .
4
4
1 0 1
− −
12 12
Thus,
2 1 1
0 − 6 0 12 − 12
(1)
1 2 1 1
e = 7 − 0 = − .
7
4
28
3 2 1 1
− 0 −
9 9 12 36
Similarly, for the second approximation of the error, we have to compute
the following:
2
0 −6 0
(2)
1 2 (1)
e = 7− 0 e
7
3 2
− 0
9 9
268 Applied Linear Algebra and Optimization using MATLAB
or
2 1 1
0 − 0 − 12 84
6
1 2 − 1 = 5 ,
e(2) =
−7 0
7 28 252
3 2 1 5
− 0
9 9 36 252
(a) Now by using the Gauss–Seidel method, first we compute the Gauss–
Seidel iteration matrix TG and the vector c as follows:
1 1
0 − 0
3
6
1 2 11
TG =
0
and c= .
21 7
42
23 4 41
0 −
189 63 378
1 1
0 − 0
3
6
1 2 (k)
11
x(k+1) =
0 x + , k = 0, 1, 2.
21 7
42
23 4 41
0 −
189 63 378
Iterative Methods for Linear Systems 269
(d) The first and second approximations of the error can be calculated as
follows:
1
0 − 0
3
1 2 e(0) = [− 1 , − 1 , 19 ]T
(1)
e = 0
21 7 12 84 756
23 4
0
189 63
and
1
0 − 0
3
1 e(1) = [ 1 , 5 , 1 ]T ,
2
e(2) =
0
21 7 252 756 6804
23 4
0
189 63
which is the required second approximation of the error. •
kx − x(k) k ≤ kT kk kx(0) − xk
(2.20)
(k) kT kk
kx − x k ≤ kx(1) − x(0) k.
1 − kT k
Note that the smaller the value of kT k, the faster the convergence of
the iterative methods.
the Gauss–Seidel iterative method converges faster than the Jacobi iterative
method.
Solution. Here we will show that the l∞ -norm of the Gauss–Seidel itera-
tion matrix TG is less than the l∞ -norm of the Jacobi iteration matrix TJ ,
i.e.,
kTG k < kTJ k.
The Jacobi iteration matrix TJ can be obtained from the given matrix A as
272 Applied Linear Algebra and Optimization using MATLAB
follows:
1
0 0
5
−1
5 0 0 0 0 −1
1
−1
TJ = −D (L+U ) = − 0 3 0 −1 0 0 = .
0 0
3
0 0 4 0 −1 0
1
0 − 0
4
Then the l∞ -norm of the matrix TJ is
1 1 1 1
kTJ k∞ = max , , = = 0.3333 < 1.
5 3 4 3
The Gauss–Seidel iteration matrix TG is defined as
5 0 0 0 0 −1
TG = −(D + L)−1 U = − −1 3 0 0 0 0 ,
0 −1 4 0 0 0
and it gives
1
0 0 5
1
TG = 0 0 −
.
15
1
0 0
60
Then the l∞ -norm of the matrix TG is
1 1 1 1
kTG k∞ = max , , = = 0.2000 < 1,
5 15 60 5
which shows that the Gauss–Seidel method will converge faster than the
Jacobi method for the given linear system. •
Note that the condition kT k < 1 is equivalent to the condition that a
matrix A is to be strictly diagonally dominant.
Iterative Methods for Linear Systems 273
For the Jacobi method for a general matrix A, the norm of the Jacobi
iteration matrix is defined as
n
X aij
kTJ k = max aii .
1≤i≤n
j=1
j6=i
TJ = −D−1 (L + U ),
1
2 1 1
0 − − −
0 0 0 10 10 10
10
0 2 1 1
1 1 2
0 1 − 0 − −
0 0
1 0 1 2 12 12 12
TJ = −
12 = ,
2 1 0 3
2 1 3
0 1
1 2 1 0
− − 0 −
0 0
13 13 13
13
1
0 0 0 1 2 1
15 − − − 0
15 15 15
then the l∞ norm of the matrix TJ is
4 4 6 4 6
kTJ k∞ = max , , , = = 0.46154 < 1.
10 12 13 15 13
Thus, the Jacobi method will converge for the given linear system.
x(1) = [0.5, 0.75, 0.07692, 0.86667]T , x(2) = [0.25564, 0.62970, −0.12436, 0.72821]T .
kTJ kk
kx − x(k) k ≤ kx(1) − x(0) k ≤ 10−4 .
1 − kTJ k
It gives
6 k
( 13 )
7 (0.86667) ≤ 10−4
13
or
7
6 k ( 13 ) × 10−4
( ) ≤ ,
13 0.86667
which gives
(0.46154)k ≤ (6.21) × 10−5 .
Taking ln on both sides, we obtain
2
k ln( ) ≤ ln (6.21) × 10−5
3
or
k(−0.77327) ≤ (−9.68621),
and it gives
k ≥ 12.5263 or k = 13,
276 Applied Linear Algebra and Optimization using MATLAB
TG = −(D + L)−1 U,
Thus, the Gauss–Seidel method will converge for the given linear system.
Iterative Methods for Linear Systems 277
kTJ kk
kx − x(k) k ≤ kx(1) − x(0) k ≤ 10−4 .
1 − kTJ k
It gives
(0.4)k
(0.74252) ≤ 10−4
0.6
or
(0.4)k ≤ (8.08 × 10−5 ).
278 Applied Linear Algebra and Optimization using MATLAB
or
k(−0.91629) ≤ (−9.4235),
and it gives
k ≥ 10.28441 or k = 11,
Example 2.10 Solve the following system of linear equations using Gauss–
Seidel iterative methods, using = 10−5 in the l∞ -norm and taking the
initial solution x(0) = [0, 0, 0, 0]T :
5x1 − x3 = 1
14x2 − x3 − x4 = 1
−x1 − x2 + 13x3 = 4
− x2 + 9x4 = 3.
(k+1) 1h (k)
i
x1 = 1 + x3
5
(k+1) 1h (k+1)
i
x4 = 3 + x2 .
9
(0) (0) (0) (0)
So starting with an initial approximation x1 = 0, x2 = 0, x3 = 0, x4 =
0, and for k = 0, we get
(1) 1h (0)
i
x1 = 1 + x3 = 0.200000
5
(1) 1h (1)
i
x4 = 3 + x2 = 0.341270.
9
The first and subsequent iterations are listed in Table 2.9. •
Note that the Gauss–Seidel method converged very fast (only five it-
erations) and the approximate solution of the given system [0.267505,
0.120302, 0.337524, 0.346700]T is equal to the exact solution [0.267505,
0.120302, 0.337524, 0.346700]T up to six decimal places.
280 Applied Linear Algebra and Optimization using MATLAB
Example 2.11 Find the eigenvalues and eigenvectors of the following ma-
trix:
−6 0 0
A = 11 −3 0 .
−3 6 7
Solution. To find the eigenvalues of the given matrix A by using (2.23),
we have
−6 − λ 0 0
11 −3 − λ 0 = 0,
−3 6 7−λ
which gives a characteristic equation of the form
λ3 + 2λ2 − 45λ − 126 = 0.
It factorizes to
(−6 − λ)(−3 − λ)(7 − λ) = 0
and gives us the eigenvalues λ = −6, λ = −3, and λ = 7 of the given
matrix A. Note that the sum of these eigenvalues is −2, and this agrees
with the trace of A. After finding the eigenvalues of the matrix we turn to
the problem of finding eigenvectors. The eigenvectors of A corresponding
to the eigenvalues λ are the nonzero vectors x that satisfy (2.22). Equiv-
alently, the eigenvectors corresponding to λ are the nonzero vectors in the
solution space of (2.22). We call this solution space the eigenspace of A
corresponding to λ.
>> D = eig(A);
>> A = [4 1 − 3; 0 0 2; 0 0 − 3];
>> B = max(eig(A))
B=
4
gives r r
cb cb
λ1 = − , λ2 =
ad ad
and r
cb
λmax = .
ad
Similarly, we can find the Gauss–Seidel iteration matrix as
TG = −(L + D)−1 U,
which gives
cb
µ1 = 0, µ2 =
ad
and
cb
µmax = .
ad
Thus,
r !2
cb
µmax = = λ2max ,
ad
which is the required result. •
286 Applied Linear Algebra and Optimization using MATLAB
The necessary and sufficient condition for the convergence of the Jacobi
iterative method and the Gauss–Seidel iterative method is defined in the
following theorem.
Note that the condition ρ(T ) < 1 is satisfied when kT k < 1 because
ρ(T ) ≤ kT k for any natural norm. •
Example 2.13 Find the spectral radius of the Jacobi and the Gauss–Seidel
iteration matrices using each of the following matrices:
2 0 −1 1 −1 1
(a) A = −1 3 0 , (b) A = −2 2 −1 ,
0 −1 4 0 1 5
1 0 0 1 0 −1
(c) A = −1 2 0 , (d) A = 1 1 0 .
0 −1 3 0 1 1
Iterative Methods for Linear Systems 287
Solution. (a) The Jacobi iteration matrix TJ for the given matrix A can
be obtained as
1
0 0
2
1
TJ =
3 0 0 ,
1
0 0
4
and the characteristic equation of the matrix TJ is
1
det(TJ − λI) = −λ3 + = 0.
24
Solving this cubic polynomial, the maximum eigenvalue (in absolute) of TJ
329
is , i.e.,
949
329
ρ(TJ ) = = 0.3467.
949
Also, the Gauss–Seidel iteration matrix TG for the given matrix A is
1
0 0
2
1
TG = 0 0
6
1
0 0
24
and has the characteristic equation of the form
1 2
det(TG − λI) = −λ3 + λ = 0.
24
Solving this cubic polynomial, we obtain the maximum eigenvalue of TG ,
1
, i.e.,
24
1
ρ(TG ) = = 0.0417.
24
288 Applied Linear Algebra and Optimization using MATLAB
and
0 0 1 0 0 1
TJ = −1 0 0 , TG = 0 0 −1 ,
0 −1 0 0 0 1
with
ρ(TJ ) = 1.0000 and ρ(TG ) = 1.0000,
respectively. •
and it gives
k
1 k
lim = 0 and lim = 0.
k→∞ 3 k→∞ 3k+1
Hence, the given matrix A is convergent. •
Since the above matrix has the eigenvalue 31 of order two, its spectral
radius is 13 . This shows the important relation existing between the spectral
radius of a matrix and the convergent of a matrix.
Theorem 2.6 The following statements are equivalent:
1. A is a convergent matrix.
2. lim kAn k = 0, for all natural norms.
n→∞
3. ρ(A) < 1.
4. lim An x = 0, for every x. •
n→∞
Solving this cubic equation, the eigenvalues of A are –2, –1, and 1. Thus
the spectral radius of A is
Also,
−2 1 0 −2 1 2 5 −2 −4
AT A = 1 0 1 1 0 0 = −2 2 2 ,
2 0 0 0 1 0 −4 2 4
which gives the eigenvalues 0.4174, 1, and 9.5826. Therefore, the spectral
radius of AT A is 9.5826. Hence,
p √
kAk2 = ρ(AT A) = 9.5826 ≈ 3.0956.
292 Applied Linear Algebra and Optimization using MATLAB
−λ3 + 4λ2 + 9λ − 36 = 0.
Note that this result is also true for any natural norm. •
which gives the eigenvalues 17.96 and 0.04. The spectral radius of AT A is
17.96. Hence, p √
kAk2 = ρ(AT A) = 17.96 ≈ 4.24.
Since a characteristic equation of (A−1 )T (A−1 ) is
−1 T −1
13 − λ 4
det[(A ) (A ) − λI] = = λ2 − 18λ + 49 = 0,
4 5−λ
294 Applied Linear Algebra and Optimization using MATLAB
which gives the eigenvalues 14.64 and 3.36 of (A−1 )T (A−1 ), its spectral
radius 14.64. Hence,
p √
kA−1 k2 = ρ((A−1 )T (A−1 )) = 14.64 ≈ 3.83.
Note that the eigenvalues of A are 3.73 and 0.27, therefore, its spectral
radius is 3.73. Hence,
1
< |3.73| < 4.24,
3.83
which satisfies Theorem 2.9. •
which is equivalent to
x(k+1) = Tω x(k) + c, (2.28)
where
Tω = (D + ωL)−1 [(1 − ω)D − ωU ] and c = ω(D − ωL)−1 b (2.29)
are called the SOR iteration matrix and the vector, respectively.
which is equal to
−1
5 0 −0.025 1.005
Tω = .
−1.005 10 0 −0.05
Thus,
0.2 0 −0.025 1.005
Tω =
0.0201 0.1 0 −0.05
or
−0.005 0.201
Tω = .
−0.0005 0.0152
The l∞ -norm of the matrix Tω is
kTω k∞ = max{0.206, 0.0157} = 0.206.
•
Example 2.20 Solve the following system of linear equations, taking an
initial approximation x(0) = [0, 0, 0, 0]T and with = 10−4 in the l∞ -norm:
2x1 + 8x2 = 1
5x1 − x2 + x3 = 2
−x1 + x2 + 4x3 + x4 = 12
x2 + x3 + 5x4 = 12.
(a) Using the Gauss–Seidel method.
(b) Using the SOR method with ω = 0.33.
(1) 1h (0)
i
x1 = 1 − 8x2 = 0.5
2
Example 2.21 Solve the following system of linear equations using the
SOR method, with = 0.5 × 10−6 in the l∞ -norm:
2x1 + x2 = 4
x1 + 2x2 + x3 = 8
x2 + 2x3 + x4 = 12
x3 + 2x4 = 11.
Start with an initial approximation x(0) = [0, 0, 0, 0]T and take ω = 1.27.
Iterative Methods for Linear Systems 299
Solution. For the given system, the SOR method with ω = 1.27 is
we obtain
(1) (0) 1.27 (0)
x1 = (1 − 1.27)x1 + [4 − x2 ] = 2.54
2
In practice, ω should be chosen in the range 1 < ω < 2, but the precise
choice of ω is a major problem. Finding the optimum value for ω depends
on the particular problem (size of the system of equations and the nature
of the equations) and often requires careful work. A detailed study for
the optimization of ω can be found in Isaacson and Keller (1966). The
following theorems can be used in certain situations for the convergence of
the SOR method.
302 Applied Linear Algebra and Optimization using MATLAB
Now to find the spectral radius of the Jacobi iteration matrix TJ , we use
the characteristic equation
λ
det(TJ − λI) = |TJ − λi| = −λ3 + ,
2
which gives the eigenvalues of matrix TJ , as λ = 0, ± √12 . Thus,
1
ρ(TJ ) = √ = 0.707107,
2
and the optimal value of ω is
2
ω= p = 1.171573.
1 + 1 − (0.707107)2
Also, note that the Gauss–Seidel iteration matrix TG has a the form
1
0 2 0
1 1
TG = 0
,
4 2
1 1
0
8 4
and its characteristic equation is
λ2
det(TG − λI) = |TG − λI| = −λ3 + .
2
Thus,
1
= 0.50000 = (ρ(TJ ))2 ,
ρ(TG ) =
2
which agrees with Theorem 2.12. •
Note that the optimal value of ω can also be found by using (2.30) if the
eigenvalues of the Jacobi iteration matrix TJ are real and 0 < ρ(TJ ) < 1. •
304 Applied Linear Algebra and Optimization using MATLAB
Example 2.23 Find the optimal choice for the relaxation factor ω by us-
ing the matrix
5 −1 −1 −1
2 5 −1 0
A= −1
.
−1 5 −1
−1 −1 −1 5
Solution. Using the given matrix A, we can find the Jacobi iteration
matrix TJ as
1 1 1
0
5 5 5
2 1
5 0 5 0
TJ = .
1 1 1
0
5 5 5
1 1 1
0
5 5 5
Now to find the spectral radius of the Jacobi iteration matrix TJ , we use
the characteristic equation
det(TJ − λI) = 0,
6 2 8 8
−λ4 − λ − λ− = (5λ − 3) ∗ (5λ + 1)3 = 0.
25 125 125
Solving the above polynomial equation, we obtain
3 1 1 1
λ = ,− ,− ,− ,
5 5 5 5
which are the eigenvalues of the matrix TJ . From this we get
3
ρ(TJ ) = = 0.6,
5
Iterative Methods for Linear Systems 305
λ4 − 0.1875λ2 + 1/256 = 0.
which shows that the Jacobi method will converge for the given linear sys-
tem.
ρ(Tω ) = ω − 1.
Now to find the spectral radius of the SOR iteration matrix Tω , we have to
calculate first the optimal value of ω by using
2
ω= p .
1 + 1 − [ρ(TJ )]2
Iterative Methods for Linear Systems 307
Using this optimal value of ω, we can compute the spectral radius of the
SOR iteration matrix Tω as follows:
Thus the SOR method will also converge for the given system, and faster
than the other two methods, because
Program 2.3
MATLAB m-file for the SOR Iterative Method
function sol=SORM(Ab,x,w,acc) % Ab = [A b]
[n,t]=size(Ab); b=Ab(1:n,t); R=1; k=1;
d(1,1:n+1)=[0 x];
k=k+1; while R > acc
for i=1:n
sum=0;
for j=1:n
if j <= i − 1; sum = sum + Ab(i, j) ∗ d(k, j + 1);
elseif j >= i + 1; sum = sum + Ab(i, j) ∗ d(k − 1, j + 1);
end;end
x(1, i) = (1 − w) ∗ d(k − 1, i + 1) + (w/Ab(i, i)) ∗ (b(i, 1) −
sum);
d(k, 1) = k − 1; d(k, i + 1) = x(1, i); end
R = max(abs((d(k, 2 : n + 1) − d(k − 1, 2 : n + 1))));
if R > 100 & k > 10; break; end
k=k+1; end; x=d;
308 Applied Linear Algebra and Optimization using MATLAB
This method has merit for nonlinear systems and optimization prob-
lems, but it is not used for linear systems because of slow convergence. An
alternative approach uses a set of nonzero direction vectors {v(1) , . . . , v(n) }
that satisfy
< v(i) , Av(j) >= 0, if i 6= j.
This is called an A-orthogonality condition, and the set of vectors {v(1) , . . . , v(n) }
is said to be A-orthogonal.
In the conjugate gradient method, we use v(1) equal to r(0) only at the
beginning of the process. For all later iterations, we choose
kr(k) k2 (k)
v(k+1) = r(k) + v
kr(k−1) k2
to be conjugate to all previous direction vectors.
Note that the initial approximation x(0) can be chosen by the user,
with x(0) = 0 as the default. The number of iterations, m ≤ n, can be
chosen by the user in advance; alternatively, one can impose a stopping
criterion based on the size of the residual vector, kr(k) k, or, alternatively,
the distance between successive iterates, kx(k+1) − x(k) k. If the process is
carried on to the bitter end, i.e., m = n, then, in the absence of round-off
errors, the results will be the exact solution to the linear system. More
iterations than n may be required in practical applications because of the
introduction of round-off errors.
Example 2.25 The linear system
2x1 − x2 = 1
−x1 + 2x2 − x3 = 0
− x2 + x3 = 1
has the exact solution x = [2, 3, 4]T . Solve the system by the conjugate
gradient method.
Solution. Start with an initial approximation x(0) = [0, 0, 0]T and find the
residual vector as
which, as one can check, is conjugate to both v(1) and v(2) . Thus, the
solution is obtained from
13 1
6 −4
2
kr (2) 2
k 1
(3) (2) (3) 4
x =x + v = 3 + 3 0 = 3 .
< v(3) , v(3) >
4
8
11 1
3 2
Since we appl