2. Now investigate products of lower triangular matrices which have all diagonal entries equal to 1. Such
a matrix is called a unit lower triangular matrix. For example you could type
n = 2
L1 = tril(rand(n),-1)+eye(n),L2 = tril(rand(n),-1)+eye(n),L1*L2,L2*L1
Execute this line several times and inspect the result each time. Repeat with n = 3,4,5.
The second input -1 for tril causes zeros to be put above the first subdiagonal, that is, on the main
diagonal as well as in the upper triangle. So adding the identity matrix eye(n) produces a random
lower triangular matrix with 1’s on the diagonal.
(a) (hand) For each pair 1
L and 2
L, what entries of 12
LL and 21
LL appear to be the same?
(b) (hand) Prove that the product of any two nn× unit lower triangular matrices will be a unit lower
triangular matrix. Notice all you need to prove here is that each diagonal entry of the product is 1 (why?).
3. (Extra Credit) Suppose L is an nn× lower triangular matrix with each diagonal entry nonzero. Create
[
]
A
LI=, where I denotes the nn× identity matrix. Explain why the reduced echelon form of A must
be of the form
[
]
IK
, where K is another nn× lower triangular matrix with nonzero diagonal entries.
MATLAB Project: The Adjacency Matrix of a Graph Name____________________________
Purpose: To learn about graph and adjacency matrix, to see how the powers of the adjacency
matrix provide information about the graph and vice versa, and to apply these ideas
Prerequisite: Section 2.1
MATLAB functions used: – , + , ^; and adjdat from Laydata4 Toolbox
Example. Verify that the graph below has the matrix A shown as its adjacency matrix.
010000
101001
⎡⎤
⎢⎥
⎢⎥
direct path from node1 to node 2
direct path from node 2 to node1; 2 to 3; 2 to 6
Theorem (Interpretation of the powers of an adjacency matrix). If
A
is the adjacency matrix of a graph,
then the (, )ijentry of k
A
is a nonnegative integer which is the number of paths of length k from node i to
node j in the graph.
1. (hand) To understand what the theorem says for the example above, let us carefully examine the (6,3)
entry of 2
A
. Using the “Row-Column Rule” from Section 2.1 in the text, the (6,3) entry of 2
A
looks like
61 13 62 23 63 33 64 43 65 53 66 63
aa aa aa aa aa aa+++++
.
Evaluate each term in this expression (you finish): (0)(0) + (1)(1) + = ________ .
5
3
2. (a) Type adjdat to get the matrix A above, then type A^2,A^3 and record the results:
2
A
= 3
A
=
(b) (hand) Note that the (1,2) entry of 2
A
is zero, so there are no paths of length two from node 1 to node 2.
Verify this by studying the graph. Similarly, notice that the (6,6) entry of 3
A
is two, so there are two paths
of length three from node 6 to itself; study the graph to see that they are 6–4–5–6 and 6–5–4–6.
In the same way, study the matrices and the graph and answer the following questions.
What are the paths of length two from node 2 to itself? ______________________________
What are the paths of length three from node 3 to node 4? ____________________________
3. (hand) Suppose A is the adjacency matrix of a graph. Explain why you must calculate the matrix sum
A
4. Eight workers, denoted W1, …, W8, handle a potentially dangerous substance. Safety precautions are
taken but accidents do happen occasionally. It is known that if a worker becomes contaminated, he or she
could spread this through contact with another worker. The following graph shows the level one contacts
between the workers.
(a) (hand) Write the adjacency matrix A for the following graph:
(b) (MATLAB) Store A, type A + A^2 + A^3 and record the result:
Use this to answer the following questions.
Which workers have contact level 3 with W3? __________________________________
Which workers have contact level 3 with W7? __________________________________
(c) Define what you mean by a worker being dangerous. Be very specific so anyone could decide
whether a worker is “dangerous” according to your definition:
(d) Which workers are the most dangerous if contaminated? ___________________
Which workers are least dangerous? _____________________
MATLAB Project: Cryptography Name_______________________________
Purpose: To see a method for using matrix operations to encode and decode messages
Prerequisite: Section 2.2 and elementary modular arithmetic
MATLAB functions used: rem, *; and crypdat from Laydata4 Toolbox
Background. Before beginning this project, read the “Notes About Arithmetic Modulo 26” at the end.
A disadvantage of a simple encoding system where each letter of the alphabet is replaced by one
other symbol is that it preserves the frequencies of individual letters. For example, the letter e is the most
common letter in English and, therefore, the letter most commonly appearing in the encoded message will
Substitution Table
ABCDEFGHIJKLMNOPQRS TUVW X YZ
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 0
1. Consider the message HURRY UP
(a) (hand) To encode this message, first divide the message into groups of two letters, adding a
“dummy” letter to fill out the last pair since the message has an odd number of letters:
HU RR YU PP
Using the Substitution Table, find the column vector for each pair of letters. The first is given, and you fill
in the rest:
(b) (MATLAB) Type crypdat to get the matrices for this project. Store the vectors above as a,
b, c and d. Transform each vector by multiplying it by the matrix 11
03
A
⎡
⎤
=
⎢
⎥
⎣
⎦, and use arithmetic modulo
26 to reduce each number you see to a new number between 0 and 25. Call the new vectors a1, b1, c1
and d1. The following line will get you started a = [8;21], A*a, a1 = rem(A*a, 26), and
⎤
⎡
8, A*a = ⎥
⎤
⎡
29 and a1 = ⎥
⎤
⎡
3.
2. (hand and MATLAB) You have received the following message from a friend:
XORLBSWOSSBUPX
Suppose that you know this message was encoded using the matrix A from question 1. To decode it you
must create a vector for each block of two letters and multiply these by a matrix B, which is the inverse of
Check the following arithmetic to verify for yourself that 117
09
B
⎡
⎤
=
⎢
⎥
⎣
⎦does satisfy these equations:
117 1 1
0903
BA ⎡⎤⎡⎤
=⎢⎥⎢⎥
⎣⎦⎣⎦
=152
027
⎡⎤
⎢⎥
⎣⎦
10
mod 26
01
⎡⎤
≡⎢⎥
⎣⎦ and 11 126 10
mod 26
03027 01
AB ⎡⎤⎡ ⎤⎡⎤
=≡
⎢⎥⎢ ⎥⎢⎥
⎣⎦⎣ ⎦⎣⎦
(a) Divide the coded message sent by your friend into blocks of two letters:
___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___
(b) Write the vector for each block:
a1 = b1 = c1 = d1 = e1 = f1 = g1 =
(c) Store your vectors a1, .., g1 and multiply each vector by B. After reducing each entry
modulo 26, record these:
a = b = c = d = e = f = g =
(d) Use the Substitution Table to retrieve the letter pairs corresponding to the vectors for the
decoded message:
a → b → c → d → e → f → g →
(e) Write the decoded message: ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___
3. (hand and MATLAB) Using the matrix A above, encode the message STAY ON THE PATH. Record
the vectors for the original message, the vectors for the coded message, and then the coded message.
5. (hand) Let 11
05
C⎡⎤
=⎢⎥
. Calculate the inverse of mod 26C. Show calculations. Call the new matrix
Remarks. The method used here is still not too safe because statistical analysis using tables of letter-pair
Page 5 of 5 MATLAB Project: Cryptography
Notes About Arithmetic Modulo 26. Two numbers r and s are called congruent modulo 26 if their
difference is an integer multiple of 26. When that is true, we write mod 26rs≡. For example,
63 11mod 26≡and 3 23mod 26−≡ are both true.
If the original number is positive, you can reduce it modulo 26 by dividing it by 26 and using the
remainder. MATLAB’s function rem does this, and it works on a single number or on a matrix. For
example, typing rem(29,26) reduces the number 29 module 26 and returns 3. If you enter x as the
A number r is said to have an inverse modulo 26 if there is another number s such that
1mod26rs ≡. For example, (3)(9) = 27 which is congruent to 1 modulo 26 so 3 and 9 are inverses
modulo 26. However, 2 has no inverse modulo 26. To see this, calculate 2 times each number between 0
and 25 and check that none of these products are congruent to 1 modulo 26. A number which has a
common factor with 26 other than 1 or -1 will not have an inverse modulo 26.
If B and C are integer matrices and reducing each number in B modulo 26 yields C, we write
mod 26BC≡. A square integer matrix C is said to have an inverse modulo 26 if there is another integer
matrix D such that mod 26CD I≡and mod 26DC I≡. It is true that such D will exist if and only if
det(C) is a number which has an inverse modulo 26. For example, the matrix A used in questions 1-4 has
determinant 3. Since the only common positive factor with 26 is 1, A has an inverse modulo 26.
To calculate the inverse of C modulo 26, calculate the reduced echelon form of
[
]
CI
by doing
10 10
2901
⎡⎤
⎢⎥
⎣⎦
~10 10
09 21
⎡⎤
⎢⎥
−
⎣⎦
~10 10
027 63
⎡⎤
⎢⎥
−
⎣⎦
10 10
mod 26
01202
⎡⎤
≡⎢⎥
⎣⎦
.
Thus, the inverse of the matrix modulo 26 is 10
20 3
⎡
⎤
⎢
⎥
⎣
⎦. Let us check to see if this is correct.
⎡
⎢
⎣
⎡
⎢
⎣
⎡
⎢
⎣
MATLAB Project: Using Backslash to Solve Ax = b Name_______________________________
Purpose: To learn about backslash and why it is the preferred method for solving systems Ax =b
when A is invertible
Prerequisite: Section 2.2 and the discussion of condition number in Section 2.3
MATLAB functions used: \, inv, format, rand, norm, diary, hilb; and ref from Laydata4
Toolbox
Background: This project is about square invertible matrices only. Suppose
A
is such a matrix and you
want to solve .A=xb
By Theorem 5 in Section 2.2, there is a unique solution to this system. MATLAB
has a special operator called backslash for solving this type of system, and it usually gives excellent
results. It is the method people use in professional settings. To use it, store A and b and type A\b .
Backslash is the best of these methods. It uses an algorithm that is fast and minimizes roundoff error. It
also checks the condition number of the coefficient matrix. If the condition number is large, it will be hard
to get an accurate answer using any numerical method. Fortunately such matrices occur rarely in real
world problems. But if backslash does detect a very large condition number, it will warn you by printing a
message “Matrix is close to singular or badly scaled. Results may be
inaccurate.” Do not ignore such a warning if you ever see it, for it means the solution is probably not
correct to very many digits. If you need more accuracy, consult a numerical analyst.
The inv function also checks the condition number, but calculating 1
A
−requires a lot more
1. (MATLAB) Here you will use the square matrices in exercises 29, 31, 39 and 41 in Section 2.2. For
each of these, you will create a linear system A=xband solve it using all three methods described above.
You will not see any warnings, so these are “good” problems. You will also see that the solutions are
almost identical as expected.
(a) To get started, determine the path where you will store your work. For example, if you install a flash
drive into the computer’s drive E: drive, type diary E:\solve to open a file called “solve” on
your flash drive.
Page 2 of 3 MATLAB Project: Using Backslash to Solve Ax = b
Type the following lines to use the matrix in exercise 29 for the problem here:
c2s2 (Opens Chapter 2 Section 2 Problems)
Type [norm(x1–x2) norm(x1–x3) norm(x2–x3)] to calculate the lengths of the differences of
the three solution vectors, and view them side by side. Record the norm of each difference vector, in the
table below. You can round each mantissa to an integer.
(b) Repeat the instructions above for the matrices in exercises 31, 39 and 41 in Section 2.2. Note that in
exercises 39 and 41 the matrix is called D, not A, so modify your commands with D\b,inv(D)*b,
rref([D b]).
Norms of the difference vectors
Exercise 29 Exercise 31 Exercise 39 Exercise 41
2. Fortunately, most matrices that show up in real world problems behave well in numerical calculations,
like those in question 1. However, it is worthwhile to see some with large condition numbers, so you
know they do exist!
One of the classic types of matrices for which none of the methods above tends to yield a very
accurate solution is the type called Hilbert matrices. There is a Hilbert matrix of every size, and
MATLAB has a special command hilb for creating them, since they are frequently used as examples
and test cases for new software.
Study these to see the pattern of entries. Think what the first and last row of the 20×20 Hilbert matrix
will look like.
3. (MATLAB) Now solve A=xb
when A is the 20×20 Hilbert matrix, using the three methods above.
You already know what this big matrix looks like, so use semicolons as indicated to avoid printing the
matrix and the large vectors:
A = hilb(20); b = rand(20,1); (create the 20×20 Hilbert matrix and a random 20×1 column)
x1 = A\b;
If you did not see warnings from the calculations A\b and inv(A)*b , you may be working on a
version which keeps more significant digits in its floating point calculations. If thisis the case , repeat the
calculations using larger sizes of Hilbert matrices until you do see warnings by typing . so try
(a) Record the size which finally produces warnings: ____________.
After solving Ax = b with a matrix that causes warnings to appear, type the following line in order
to view the three solutions you have created, in a format that makes it fairly easy to compare them:
(b) The mantissas of the entries in x1 and x2 may agree in some digits, but this is quite misleading.
To see how different these vectors really are, type [norm(x1–x2 norm(x1–x3) norm(x2–x3)]
to calculate the norms of the difference vectors, and record those. Round to integers as before:
Norms of the difference vectors for the Hilbert matrix problem
x1–x2: ________________ x1–x3: ________________ x2–x3: ________________
Clearly the warnings are justified! Even the difference x1-x2 is very large, and x1 and x2 were
4. To finish, type diary off, which will close the file called “solve” on your computer. Exit
MATLAB and open this file with your favorite text editor. If it contains more than 4 pages, try to reduce
its length before printing. For example, erase unnecessary blank lines or big matrices you may have
created, and perhaps reduce the font size. Print the file and attach that printout to this project.
Caution: This project is about square invertible coefficient matrices only. Do not use backslash
MATLAB Project: Roundoff Error in Matrix Calculations Name______________________
Purpose: To see some of the effects of roundoff error by solving a variety of linear systems, using
three different algorithms
Prerequisite: Section 2.2 and the discussion of condition number in Section 2.3
MATLAB functions used: *, –, \, format, diary, inv, rand, norm, hilb, cond; and ref from
Laydata4 Toolbox
Background. Most computer calculations use floating point arithmetic, not exact arithmetic. This means
each number is stored in “mantissa-exponent” format, with only the first few significant digits kept in the
mantissa and all trailing digits simply lost, or “chopped.” In base 10 a floating point number looks like
Computers usually use base 2 arithmetic but display numbers in base 10 format. Also, they have
more digits stored for a number than they usually display. You can decide how many digits you want to
see, which exponents are shown, or other displays by changing the display format. Changing the display
format does not change the number that is stored, only the way it gets printed out.
1. (MATLAB) Type the following commands to see some of the different ways MATLAB can display
numbers. It starts out in the default mode, which is called “format short.” For each command, record the
result on the right:
x = 12.123456789123456789
format long, x
y = .0012123456789123456789
format long, y
Page 2 of 6 MATLAB Project: Roundoff Error in Matrix Calculations
More background. It is important to be able to estimate the accuracy of a computer solution and to
detect when a solution is likely to be inaccurate. Here is one way to check accuracy. When 1
xis a
calculated solution to A=xb
, the vector 11
A=−rb x
is called the residual vector, or just the residual.
Theoretically, 1
r should be identically zero, but usually it is not because of roundoff error. However, if 1
r
is quite small, that is a good sign that the calculated solution 1
xis probably reasonably accurate.
It is helpful to define the norm, or length of a vector x as 22 2
12 n
x
xx=+++x…. For example, if
[
]
42 21=− −x, then 5=x.
A remarkable fact about professionally written matrix software today is that the algorithms used do
give very accurate answers for the vast majority of real world problems. The algorithms in such software
are also efficient, meaning they do a minimum amount of calculation and run relatively fast, even on very
big problems.
However, there are some matrices for which computer calculations are likely to be inaccurate, no
matter how good the algorithm. Professional software checks and when it detects such a matrix, it warns
the user that some matrix involved in the calculation is poorly conditioned. For now, think of
A
being
poorly conditioned as meaning the entries of
A
and 1
A
−differ dramatically in size. The opposite of being
poorly conditioned is being well conditioned. There is discussion about condition number in Section 2.3
One more definition will be useful. A Hilbert matrix is of the form
112 1
12 13 1( 1)
⎢
⎣
n
n
⎡
⎤
⎢
⎥
+
⎢
⎥
⎢
⎥
.
Page 3 of 6 MATLAB Project: Roundoff Error in Matrix Calculations
Instructions: To begin, type format short, format compact . That second command will
suppress some of the blank lines that usually appear between MATLAB calculations allowing more
information to fit on the screen at the same time.
2. (MATLAB) Type the following lines to create 8 matrices and 4 vectors. The command rand(n,m)
creates an nm×matrix with random number entries. For the larger matrices, use semicolons as indicated
so you do not have to look at so much data:
3. (MATLAB) Solve the equation 11
A=xb
, using three different methods as indicated below. Each
method is theoretically correct but you will see that the solutions differ, especially for the larger size
matrices. The three methods are discussed on page 7 below.
(a) For simplicity, rename the matrix and vector and store the size by typing
Then type the following three lines to solve the equation A=xb
three different ways and to calculate the
residual vector for each method. Notice semicolons are used to suppress printing the vectors at first. Then
the solutions are printed side by side followed by the residual vectors printed side by side. For the larger
vectors you will appreciate the economy of space and how easy it is to compare them when displayed this
way.
x1 = A\b; res1 = b – A*x1;
(b) Type the following line and inspect to find the power of 10 that appears in these norms. In the
A1 column of Table 1 on the next page, record these powers of 10.
(c) Repeat the calculations above for 11
H=xb
. First type A = H1 into MATLAB. To avoid
having to retype all the other MATLAB commands, just press the up arrow key until you find the line you
want next, and press [Enter] to execute it again. Record your powers of 10 in Table 1 on the next page.
4. (MATLAB) Repeat the instructions in question 3 for 22
A=xb
and 22
H=xb
; 33
A=xb
and 33
H=xb
;
and 44
A=xb
and 44
H=xb
. For example, 2
A
has a different shape than 1
A
so you need to first type
and then do the calculations. Again you can use the up arrow key to avoid having to retype the other
commands.
Table 1. Power of 10 in the norm of each residual vector
A1 A2 A3 A4 H1 H2 H3 H4
Table 2. Warnings about poor conditioning?
A1 A2 A3 A4 H1 H2 H3 H4
x1
x2
This completes the numerical calculations for this project. To answer the rest of the questions, you will
need to inspect your diary file.
6. Examine the residual vectors res1 and res2 for each coefficient matrix. Which coefficient
matrices yielded smaller residual vectors, and which ones yielded larger residuals? Give your two lists:
7. Discuss the matrices you listed in questions 5 and 6. Are there some matrices that you think should
have triggered warnings about poor conditioning but did not? Which ones were they and why do you
think that?
Notes about the three solution methods. The following notes compare the three methods used above.
Notice we are discussing only the solution of systems Ax = b when A is square and invertible.
Page 6 of 6 MATLAB Project: Roundoff Error in Matrix Calculations
Notes about Hilbert matrices . The nn×Hilbert matrix is
112 1
12 13 1( 1)
11(1) 1/(21)
n
n
H
nn n
⎡
⎤
⎢
⎥
+
⎢
⎥
=
⎢
⎥
⎢
⎥
+−
⎣
⎦
.
Theoretically, every Hilbert matrix is invertible, but they are notoriously bad for floating point
calculations. The bigger Hilbert matrix H you use, the less accurate will be the solution to H=xb
when
calculated with floating point arithmetic. Some algorithms may appear to give somewhat more accurate
results than others, but none will be very good for a big Hilbert matrix. The problem is in the matrix
itself. Hilbert matrices are classic examples of poorly conditioned matrices which is why people use them
as test cases for algorithms and why MATLAB has the command hilb for generating them.
MATLAB actually has a function invhilb for calculating an exact inverse for a Hilbert matrix!
You might enjoy comparing invhilb(n) and inv(hilb(n)) for several values of n. The inv
function does floating point calculations, so it does not produce an exact inverse. In fact, when n is around
12 or larger, you will see that inv(hilb(n)) produces a matrix which has very few entries even close
to what they should be. Type help hilb and help invhilb for more information.
A large condition number implies that small errors in b can cause large errors in x. Roughly, if the
entries of A and b are accurate to about r significant digits and if the condition number of A is
approximately 10k, then the computed solution of A=xb
should usually be accurate to at least rk−
significant digits. A large condition number has the potential effect of losing k significant digits rendering
the solution meaningless!
MATLAB Project: Partitioned Matrices Name_______________________________
Purpose: To investigate multiplication of partitioned matrices
Prerequisite: Section 2.4
MATLAB functions used: + , – , * , eye, inv; and partdat from Laydata4 Toolbox
1. (a) (MATLAB) To get the matrices below, type partdat. (They will be named A11, A12, etc.)
⎢
⎣
⎢
⎢
⎢
⎣
⎢
⎢
⎢
⎣
⎡
⎡
⎢
12
41 1
⎡⎤
2
⎡
⎤
20
⎡
⎤
11 2
⎡
⎤
1
⎡⎤
12
⎢
⎥
⎣
⎦, 21
42
⎢
⎥
⎢
⎥
−
⎣
⎦
⎢
⎢
⎢
⎣
11 12 13
AAA
⎡⎤
11
C
⎡
⎤
(b) Calculate AC using partitioned multiplication. The following lines calculate each block.
A11*C11 + A12*C21 + A13*C31 Result =
Now write the matrix AC and mark the partitions corresponding to the blocks of A and C.