Lecture 8

Matrices and Matrix Operations in Matlab

Matrix operations

Recall how to multiply a matrix A times a vector v:

� �� � � � � �
1 2 −1 1 · (−1) + 2 · 2 3
Av = = = .
3 4 2 3 · (−1) + 4 · 2 5

This is a special case of matrix multiplication. To multiply two matrices, A and B you proceed as follows:
� �� � � � � �
1 2 −1 −2 −1 + 4 −2 + 2 3 0
AB = = = .
3 4 2 1 −3 + 8 −6 + 4 5 −2

Here both A and B are 2 × 2 matrices. Matrices can be multiplied together in this way provided that the
number of columns of A match the number of rows of B. We always list the size of a matrix by rows, then
columns, so a 3 × 5 matrix would have 3 rows and 5 columns. So, if A is m × n and B is p × q, then we can
multiply AB if and only if n = p. A column vector can be thought of as a p × 1 matrix and a row vector as
a 1 × q matrix. Unless otherwise specified we will assume a vector v to be a column vector and so Av makes
sense as long as the number of columns of A matches the number of entries in v.

Printing matrices on the screen takes up a lot of space, so you may want to use
format compact

Enter a matrix into Matlab with the following syntax:

A = [ 1 3 -2 5 ; -1 -1 5 4 ; 0 1 -9 0]

Also enter a vector u:

u = [ 1 2 3 4] ’

To multiply a matrix times a vector Au use *:

A*u

Since A is 3 by 4 and u is 4 by 1 this multiplication is valid and the result is a 3 by 1 vector.

Now enter another matrix B using

B = [3 2 1; 7 6 5; 4 3 2]

34
 35

You can multiply B times A with

B*A

but A times B is not defined and

A*B

will result in an error message.

You can multiply a matrix by a scalar:

2* A

Adding matrices A + A will give the same result:

A + A

You can even add a number to a matrix:

A + 3 % add 3 to every entry of A

Component-wise operations

Just as for vectors, adding a ’.’ before ‘*’, ‘/’, or ‘^’ produces entry-wise multiplication, division and
exponentiation. If you enter
B*B

the result will be BB, i.e. matrix multiplication of B times itself. But, if you enter
B .* B

the entries of the resulting matrix will contain the squares of the same entries of B. Similarly if you want B
multiplied by itself 3 times then enter
B ^3

but, if you want to cube all the entries of B then enter

B .^3

Note that B*B and B^3 only make sense if B is square, but B.*B and B.^3 make sense for any size matrix.

The identity matrix and the inverse of a matrix

The n × n identity matrix is a square matrix with ones on the diagonal and zeros everywhere else. It is called
the identity because it plays the same role that 1 plays in multiplication, i.e.

AI = A, IA = A, Iv = v
36 LECTURE 8. MATRICES AND MATRIX OPERATIONS IN MATLAB

for any matrix A or vector v where the sizes match. An identity matrix in Matlab is produced by the
command
I = eye (3)

A square matrix A can have an inverse which is denoted by A−1 . The definition of the inverse is that

AA−1 = I and A−1 A = I.

In theory an inverse is very important, because if you have an equation

Ax = b

where A and b are known and x is unknown (as we will see, such problems are very common and important)
then the theoretical solution is
x = A−1 b.
We will see later that this is not a practical way to solve an equation, and A−1 is only important for the
purpose of derivations. In Matlab we can calculate a matrix’s inverse very conveniently:
C = randn (5 ,5)
inv ( C )

However, not all square matrices have inverses:

D = ones (5 ,5)
inv ( D )

The “Norm” of a matrix

For a vector, the “norm” means the same thing as the length (geometrically, not the number of entries).
Another way to think of it is how far the vector is from being the zero vector. We want to measure a matrix
in much the same way and the norm is such a quantity. The usual definition of the norm of a matrix is

Definition 1 Suppose A is a m × n matrix. The norm of A is

�A� ≡ max �Av�.

�v�=1

The maximum in the definition is taken over all vectors with length 1 (unit vectors), so the definition means
the largest factor that the matrix stretches (or shrinks) a unit vector. This definition seems cumbersome at
first, but it turns out to be the best one. For example, with this definition we have the following inequality
for any vector v:
�Av� ≤ �A��v�.
In Matlab the norm of a matrix is obtained by the command
norm ( A )

For instance the norm of an identity matrix is 1:

 37

norm ( eye (100))

and the norm of a zero matrix is 0:

norm ( zeros (50 ,50))

For a matrix the norm defined above and calculated by Matlab is not the square root of the sum of the
square of its entries. That quantity is called the Froebenius norm, which is also sometimes useful, but we
will not need it.

Some other useful commands

C = rand(5,5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . random matrix with uniform distribution in [0, 1].

size(C) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . gives the dimensions (m × n) of C.
det(C) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . the determinant of the matrix.
max(C) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . the maximum of each column.
min(C) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . the minimum in each column.
sum(C) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . sums each column.
mean(C) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .the average of each column.
diag(C) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . just the diagonal elements.
C’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . tranpose the matrix.

In addition to ones, eye, zeros, rand and randn, Matlab has several other commands that automatically
produce special matrices:
hilb(6)
pascal(5)

Exercises
8.1 Enter the matrix M by
M = [1 ,3 , -1 ,6;2 ,4 ,0 , -1;0 , -2 ,3 , -1; -1 ,2 , -5 ,1]

and also the matrix  

−1 −3 3
 2 −1 6 

N = .
1 4 −1 
2 −1 2
Multiply M and N using M * N. Can the order of multiplication be switched? Why or why not? Try
it to see how Matlab reacts.
8.2 By hand, calculate Av, AB, and BA for:
    
2 4 −1 0 −1 −1 3
A =  −2 1 9 , B =  1 0 2 , v =  1 .
−1 −1 0 −1 −2 0 −1

Check the results using Matlab. Think about how fast computers are. Turn in your hand work.
8.3 (a) Write a well-commented Matlab function program myinvcheck that
38 LECTURE 8. MATRICES AND MATRIX OPERATIONS IN MATLAB

• makes a n × n random matrix (normally distributed, A = randn(n,n)),

• calculates its inverse (B = inv(A)),
• multiplies the two back together,
• calculates the residual (difference from the desired n × n identity matrix eye(n)), and
• returns the norm of the residual.
(b) Write a well-commented Matlab script program that calls myinvcheck for n = 10, 20, 40, . . . , 2i 10
for some moderate i, records the results of each trial, and plots the error versus n using a log
plot. (See help loglog.)
What happens to error as n gets big? Turn in a printout of the programs, the plot, and a very brief
report on the results of your experiments.
Lecture 9
Introduction to Linear Systems

How linear systems occur

Linear systems of equations naturally occur in many places in engineering, such as structural analysis,
dynamics and electric circuits. Computers have made it possible to quickly and accurately solve larger
and larger systems of equations. Not only has this allowed engineers to handle more and more complex
problems where linear systems naturally occur, but has also prompted engineers to use linear systems to
solve problems where they do not naturally occur such as thermodynamics, internal stress-strain analysis,
fluids and chemical processes. It has become standard practice in many areas to analyze a problem by
transforming it into a linear systems of equations and then solving those equation by computer. In this way,
computers have made linear systems of equations the most frequently used tool in modern engineering.

In Figure 9.1 we show a truss with equilateral triangles. In Statics you may use the “method of joints” to
each node of the truss1 . This set of equations is an example of a linear system. Making
write equations for √
the approximation 3/2 ≈ .8660, the equations for this truss are

.5 T1 + T2 = R1 = f1
.866 T1 = −R2 = −.433 f1 − .5 f2
−.5 T1 + .5 T3 + T4 = −f1
.866 T1 + .866 T3 = 0 (9.1)
−T2 − .5 T3 + .5 T5 + T6 = 0
.866 T3 + .866 T5 = f2
−T4 − .5 T5 + .5 T7 = 0,

where Ti represents the tension in the i-th member of the truss.

You could solve this system by hand with a little time and patience; systematically eliminating variables and
substituting. Obviously, it would be a lot better to put the equations on a computer and let the computer
solve it. In the next few lectures we will learn how to use a computer effectively to solve linear systems.
The first key to dealing with linear systems is to realize that they are equivalent to matrices, which contain
numbers, not variables.

As we discuss various aspects of matrices, we wish to keep in mind that the matrices that come up in
engineering systems are really large. It is not unusual in real engineering to use matrices whose dimensions
are in the thousands! It is frequently the case that a method that is fine for a 2 × 2 or 3 × 3 matrix is totally
inappropriate for a 2000 × 2000 matrix. We thus want to emphasize methods that work for large matrices.
1
See http://en.wikipedia.org/wiki/Truss or http://en.wikibooks.org/wiki/Statics for reference.

39
40 LECTURE 9. INTRODUCTION TO LINEAR SYSTEMS

f1 ✲ � �
4
B D

1 3 5 7

✉A � E ✉
C
2 6
R1 ✛
✻ ❄ ✻
R2 f2 R3

Figure 9.1: An equilateral truss. Joints or nodes are labeled alphabetically, A, B, . . . and Members
(edges) are labeled numerically: 1, 2, . . . . The forces f1 and f2 are applied loads and R1 , R2 and
R3 are reaction forces applied by the supports.

Linear systems are equivalent to matrix equations

The system of linear equations

x1 − 2x2 + 3x3 = 4
2x1 − 5x2 + 12x3 = 15
2x2 − 10x3 = −10
is equivalent to the matrix equation
    
1 −2 3 x1 4
 2 −5 12   x2  =  15  ,
0 2 −10 x3 −10
which is equivalent to the augmented matrix
 
1 −2 3 4
 2 −5 12 15  .
0 2 −10 −10
The advantage of the augmented matrix, is that it contains only numbers, not variables. The reason this is
better is because computers are much better in dealing with numbers than variables. To solve this system,
the main steps are called Gaussian elimination and back substitution.

The augmented matrix for the equilateral truss equations (9.1) is given by
 
.5 1 0 0 0 0 0 f1
 .866 0 0 0 0 0 0 −.433f 1 − .5 f2

 
 −.5 0 .5 1 0 0 0 −f 
 1 
 .866 0 .866 0 0 0 0 0 . (9.2)
 
 0 −1 −.5 0 .5 1 0 0 
 
 0 0 .866 0 .866 0 0 f2 
0 0 0 −1 −.5 0 .5 0
Notice that a lot of the entries are 0. Matrices like this, called sparse, are common in applications and there
are methods specifically designed to efficiently handle sparse matrices.
 41

Triangular matrices and back substitution

Consider a linear system whose augmented matrix happens to be

 
1 −2 3 4
 0 −1 6 7  . (9.3)
0 0 2 4

Recall that each row represents an equation and each column a variable. The last row represents the equation
2x3 = 4. The equation is easily solved, i.e. x3 = 2. The second row represents the equation −x2 + 6x3 = 7,
but since we know x3 = 2, this simplifies to: −x2 + 12 = 7. This is easily solved, giving x2 = 5. Finally,
since we know x2 and x3 , the first row simplifies to: x1 − 10 + 6 = 4. Thus we have x1 = 8 and so we know
the whole solution vector: x = �8, 5, 2�. The process we just did is called back substitution, which is both
efficient and easily programmed. The property that made it possible to solve the system so easily is that
A in this case is upper triangular. In the next section we show an efficient way to transform an augmented
matrix into an upper triangular matrix.

Gaussian Elimination

Consider the matrix  

1 −2 3 4
A=  2 −5 12 15  .
0 2 −10 −10
The first step of Gaussian elimination is to get rid of the 2 in the (2,1) position by subtracting 2 times the
first row from the second row, i.e. (new 2nd = old 2nd - (2) 1st). We can do this because it is essentially the
same as adding equations, which is a valid algebraic operation. This leads to
 
1 −2 3 4
 0 −1 6 7 .
0 2 −10 −10

There is already a zero in the lower left corner, so we don’t need to eliminate anything there. To eliminate
the third row, second column, we need to subtract −2 times the second row from the third row, (new 3rd =
old 3rd - (-2) 2nd), to obtain  
1 −2 3 4
 0 −1 6 7  .
0 0 2 4
This is now just exactly the matrix in equation (9.3), which we can now solve by back substitution.

Matlab’s matrix solve command

In Matlab the standard way to solve a system Ax = b is by the command

x = A \ b

This command carries out Gaussian elimination and back substitution. We can do the above computations
as follows:
42 LECTURE 9. INTRODUCTION TO LINEAR SYSTEMS

A = [1 -2 3 ; 2 -5 12 ; 0 2 -10]
b = [4 15 -10] ’
x = A \ b

Next, use the Matlab commands above to solve Ax = b when the augmented matrix for the system is
 
1 2 3 4
 5 6 7 8 ,
9 10 11 12

by entering
x1 = A \ b

Check the result by entering

A * x1 - b

You will see that the resulting answer satisfies the equation exactly. Next try solving using the inverse of A:
x2 = inv ( A )* b

This answer can be seen to be inaccurate by checking

A * x2 - b

Thus we see one of the reasons why the inverse is never used for actual computations, only for theory.

Exercises

9.1 Set f1 = 1000N and f2 = 5000N in the equations (9.1) for the equailateral truss. Input the coefficient
matrix A and the right hand side vector b in (9.2) into Matlab. Solve the system using the command
\ to find the tension in each member of the truss. Save the matrix A as A_equil_truss and keep it
for later use. (Enter save A_equil_truss A.) Print out and turn in A, b and the solution x.
9.2 Write each system of equations as an augmented matrix, then find the solutions using Gaussian
elimination and back substitution (the algorithm in this chapter). Check your solutions using Matlab.

(a)

x1 + x2 = 2
4x1 + 5x2 = 10

(b)

x1 + 2x2 + 3x3 = −1
4x1 + 7x2 + 14x3 = 3
x1 + 4x2 + 4x3 = 1
Lecture 10
Some Facts About Linear Systems

Some inconvenient truths

In the last lecture we learned how to solve a linear system using Matlab. Input the following:
A = ones (4 ,4)
b = randn (4 ,1)
x = A \ b

As you will find, there is no solution to the equation Ax = b. This unfortunate circumstance is mostly the
fault of the matrix, A, which is too simple, its columns (and rows) are all the same. Now try
b = ones (4 ,1)
x = [ 1 0 0 0] ’
A*x

So the system Ax = b does have a solution. Still unfortunately, that is not the only solution. Try
x = [ 0 1 0 0] ’
A*x

We see that this x is also a solution. Next try

x = [ -4 5 2.27 -2.27] ’
A*x

This x is a solution! It is not hard to see that there are endless possibilities for solutions of this equation.

Basic theory

The most basic theoretical fact about linear systems is

Theorem 1 A linear system Ax = b may have 0, 1, or infinitely many solutions.

Obviously, in most engineering applications we would want to have exactly one solution. The following two
theorems tell us exactly when we can and cannot expect this.

43
44 LECTURE 10. SOME FACTS ABOUT LINEAR SYSTEMS

Theorem 2 Suppose A is a square (n × n) matrix. The following are all equivalent:

1. The equation Ax = b has exactly one solution for any b.

2. det(A) �= 0.
3. A has an inverse.
4. The only solution of Ax = 0 is x = 0.
5. The columns of A are linearly independent (as vectors).
6. The rows of A are linearly independent.

If A has these properties then it is called non-singular.

On the other hand, a matrix that does not have these properties is called singular.

Theorem 3 Suppose A is a square matrix. The following are all equivalent:

1. The equation Ax = b has 0 or ∞ many solutions depending on b.

2. det(A) = 0.
3. A does not have an inverse.
4. The equation Ax = 0 has solutions other than x = 0.
5. The columns of A are linearly dependent as vectors.
6. The rows of A are linearly dependent.

To see how the two theorems work, define two matrices (type in A1 then scroll up and modify to make A2) :
   
1 2 3 1 2 3
A1 =  4 5 6 , A2 =  4 5 6 ,
7 8 9 7 8 8

and two vectors:    

0 1
b1 =  3  , b2 =  3  .
6 6
First calculate the determinants of the matrices:
det ( A1 )
det ( A2 )

Then attempt to find the inverses:

inv ( A1 )
inv ( A2 )

Which matrix is singular and which is non-singular? Finally, attempt to solve all the possible equations
Ax = b:
 45

x = A1 \ b1
x = A1 \ b2
x = A2 \ b1
x = A2 \ b2

As you can see, equations involving the non-singular matrix have one and only one solution, but equation
involving a singular matrix are more complicated.

The Residual

Recall that the residual for an approximate solution x of an equation f (x) = 0 is defined as r = �f (x)�. It
is a measure of how close the equation is to being satisfied. For a linear system of equations we define the
residual vector of an approximate solution x by

r = Ax − b.

If the solution vector x were exactly correct, then r would be exactly the zero vector. The size (norm) of r
is an indication of how close we have come to solving Ax = b. We will refer to this number as the scalar
residual or just Residual of the approximation solution:

r = �Ax − b�. (10.1)

Exercises

10.1 By hand, find all the solutions (if any) of the following linear system using the augmented matrix and
Gaussian elimination:
x1 + 2x2 + 3x3 = 4,
4x1 + 5x2 + 6x3 = 10,
7x1 + 8x2 + 9x3 = 14 .

10.2 (a) Write a well-commented Matlab function program mysolvecheck with input a number n that
makes a random n × n matrix A and a random vector b, solves the linear system Ax = b,
calculates the scalar residual r = �Ax − b�, and outputs that number as r.
(b) Write a well-commented Matlab script program that calls mysolvecheck 10 times each for
n = 10, 50, 100, 500, and 1000, records and averages the results and makes a log-log plot of the
average e vs. n. Increase the maximum n until the program stops running within 5 minutes.
Turn in the plot and the two programs.