1

MECE 390-A1
Numerical Methods of Mechanical Engineering

Chapter 3. Systems of linear equations

Lecture 7. 2018/09/25

Department of Mechanical Engineering, University of Alberta

MEC E 390 -A1

Review 2

A vector 𝒂 can be written in the rectangular cartesian coordinate:

𝒂 = 𝑎𝑥 𝒊 + 𝑎𝑦 𝒋 + 𝑎𝑧 𝒌.
Here, we refer 𝑎𝑥 , 𝑎𝑦 , 𝑎𝑧 as components of a row vector 𝒂. 𝑎𝑥 , 𝑎𝑦 , 𝑎𝑧 is also considered as [1 × 3]
matrix. Similarly, we could think of a column vector as:
𝑏𝑥 𝒊
𝒃 = 𝑏𝑦 𝒋
𝑏𝑧 𝒌
The components of the above is given by
𝑏𝑥
𝑏𝑦 ⇒ 3 × 1 matrix
𝑏𝑧
Lastly, we could think of higher order vector also known as ‘tensor’ (e.g. stress tensor). The
components of a tensor is referred as matrix. For example, components of a stress tensor can be
𝜎11 𝜎12 𝜎13
expressed as 𝜎𝑖𝑗 = 𝜎21 𝜎22 𝜎23 3 × 3 matrix.
𝜎31 𝜎32 𝜎33

MEC E 390 -A1

Review 3

In this course, we denote components of a row vector, column vector and tensor respectively as
𝑏1 𝑎11 𝑎12 𝑎1𝑛
𝑎21 𝑎22 ⋯ 𝑎2𝑛
𝑏
𝒂 = [𝑎𝑖 ] = 𝑎1 , 𝑎2 , … 𝑎𝑛 , 𝒃 = [𝑏𝑖 ] = 2 and 𝒂 = [𝑎𝑖𝑗 ] = ⋮ ⋮ ⋱ ⋮
⋮
𝑏𝑛 𝑎𝑛1 𝑎𝑛2 ⋯ 𝑎𝑛𝑛
 Matrix operations
a. Addition: 𝑨 = 𝑩 + 𝑪, 𝑨 + 𝑩 = 𝑩 + 𝑨: commutative, 𝑨 + 𝑩 + 𝑪 = 𝑨 + 𝑩 + 𝑪 : Associative.
b. Scalar multiplication: 𝑘𝑨 = 𝑩 ⇒ 𝑘𝐴𝑖𝑗 = 𝐵𝑖𝑗 (multiplying a scalar 𝑘 on each component of a
matrix i.e. 𝑘𝐴12 = 𝐵12 etc…
c. Matrix multiplication 𝑨 = 𝑩𝑪 ⇒ 𝑎𝑖𝑗 = 𝑏𝑖𝑘 𝑐𝑘𝑗 (summation over 𝑘 for 𝑘 = 1,2, … 𝑛).

Note: Matrix multiplication is NOT commutative (i.e. 𝑨𝑩 ≠ 𝑩𝑪)

There is other types of matrix multiplications e.g. Hadamard multiplication (i.e. 𝑎𝑖𝑗 = 𝑏𝑖𝑗 𝑐𝑖𝑗 ).
In this course (and in general) matrix multiplication refers to 𝑎𝑖𝑗 = 𝑏𝑖𝑘 𝑐𝑘𝑗 unless otherwise specified.

MEC E 390 -A1

Review 4

An identity matrix is defined by

1 0 ⋯ 0
0 1 ⋯ 0 1 𝑖=𝑗
𝑰 = [𝛿𝑖𝑗 ] = ⇒ 𝛿𝑖𝑗 =
⋮ ⋮ ⋱ ⋮ 0 𝑖≠𝑗
0 0 ⋯ 1
An identity matrix also satisfies
𝑨𝑰 = 𝑨.
Using the properties of an identity matrix, we define an inverse matrix as
𝑨𝑩 = 𝑰 ⇒ 𝑨−𝟏 𝑨𝑩 = 𝑨−𝟏 𝑰 ⇒ 𝑨−𝟏 𝑨 𝑩 = 𝑨−𝟏 ⇒ 𝑰𝑩 = 𝑨−𝟏 ⇒ 𝑩 = 𝑨−𝟏
The transpose of a matrix is defined as
𝑇
𝐴𝑖𝑗 = 𝐴𝑗𝑖
which can be obtained by exchanged upper and lower diagonal components of a matrix. In
particular, for a symmetric matrix, we find
𝐴𝑖𝑗 = 𝐴𝑗𝑖

MEC E 390 -A1

Systems of linear algebraic equations 5

 Many complex engineering problems can be expressed as set of equations with many
variables. Solving all equations simultaneously is equivalent to finding the roots of all equations
with respect to all variables.
To study the numerical solution of systems of linear equations, let us define more precisely our
matrix notation. Given a system of 𝒏 equation and 𝒏 unknowns:

We can write this system using matrix notation as

𝐴 𝑥 = 𝑏 , where
𝑎11 𝑎12 𝑎1𝑛 𝑥1 𝑏1
𝑎21 𝑎22 ⋯ 𝑎 𝑥2
2𝑛 𝑏
𝐴= ⋮ ⋮ ⋱ ⋮ , 𝑥 = ⋮ , 𝑏= 2
⋮
𝑎𝑛1 𝑎𝑛2 ⋯ 𝑎𝑛𝑛 𝑥𝑛 𝑏𝑛

MEC E 390 -A1

Special matrixes 6

 Examples of matrixes

𝑑11 0 ⋯ 0 a. diagonal matrix

0 𝑑22 ⋯ 0 e.g. Identity matrix is a special case of a diagonal matrix where
𝐷=
⋮ ⋮ ⋱ ⋮ 𝑑11 = 𝑑22 = ⋯ = 𝑑𝑛𝑛 = 1.
0 0 ⋯ 𝑑𝑛𝑛

b. Lower triangular 𝑢11 𝑢12 ⋯ 𝑢𝑛1 c. Upper triangular

𝑙11 0 ⋯ 0
0 𝑢22 ⋯ 0
𝑙 𝑙22 ⋯ 0 matrix 𝑈= ⋮ matrix
𝐿 = 21 ⋮ ⋱ ⋮
⋮ ⋮ ⋱ ⋮
0 0 ⋯ 𝑢𝑛𝑛
𝑙𝑛1 𝑙𝑛2 ⋯ 𝑙𝑛𝑛

𝑡11 𝑡12 0 0 ⋯ 0 d. Triangular matrix (non-zero elements appear only

𝑡21 𝑡22 𝑡23 0 ⋯ 0 along, just above or just below main diagonal.
𝑇= 0 0 ⋱ 0 0 0
⋮ 0 𝑡𝑛−1,𝑛−2 𝑡𝑛−1,𝑛−2 𝑡𝑛−1,𝑛
0 0 ⋯ 0 𝑡𝑛,𝑛−1 𝑡𝑛,𝑛

MEC E 390 -A1

Cramer’s Rule 7

Now go back to the generic matrix equation 𝐴𝒙 = 𝒃. How do we solve for 𝒙 ?

Frist of all, there is no such thing as “division by [𝐴].”

 Many different possibilities exist including:

1. Cramer’s Rule
2. Calculation of 𝐴−1
3. Gaussian Elimination
4. Iterative schemes etc…

 Cramer’s Rule
𝑎11 𝑎12 𝑥1 𝑏1 𝑎11 𝑥1 + 𝑎12 𝑥2 = 𝑏1
Consider the 2 × 2 systems of equations 𝑎 𝑎22 𝑥2 = ⇒
21 𝑏2 𝑎21 𝑥1 + 𝑎22 𝑥2 = 𝑏2
𝑏2 −𝑎21 𝑥1
Eliminating 𝑥1 from latter equations shows that 𝑥2 = .
𝑎22

MEC E 390 -A1

Cramer’s Rule 8
𝑎
Plugging this results into former equation yields 𝑎11 𝑥1 + 𝑎12 𝑏2 − 𝑎21 𝑥1 = 𝑏1 . Solving for 𝑥1 gives
22
𝑎11 𝑎22 − 𝑎12 𝑎21 𝑥1 = 𝑎22 𝑏1 − 𝑎12 𝑏2 , where
𝑎11 𝑎22 − 𝑎12 𝑎21 = determinant of 𝐴 (det(𝐴) or 𝑑𝑒𝑡𝐴 or 𝐴 for short).
𝑏1 𝑎12
Notice also that RHS = 𝑎22 𝑏2 − 𝑎12 𝑏2 = .
𝑏2 𝑎22

Thus 𝑑𝑒𝑡𝐴 𝑥1 = 𝑑𝑒𝑡𝐴1 , where 𝐴1 is obtained from 𝐴 by replacing the first column with the vector 𝒃.
𝑏 𝑎12 𝑎 𝑏1
i.e. 𝐴1 = 1 . Similarly, 𝑑𝑒𝑡𝐴 𝑥2 = 𝑑𝑒𝑡𝐴2 where 𝐴2 = 11 .
𝑏2 𝑎22 𝑎21 𝑏2

Generalizing to an 𝑛 × 𝑛 systems of linear equations gives for 1 ≤ 𝑘 ≤ 𝑛

𝑑𝑒𝑡𝐴𝑘
⇒ 𝑥𝑘 =
𝑑𝑒𝑡𝐴

Note: in the above, the solution (𝑥𝑘 ) exist only if 𝑑𝑒𝑡𝐴 ≠ 0. For 𝑑𝑒𝑡𝐴 = 0 ⇒ Singularity (no solution!!)

MEC E 390 -A1

Determinant of a matrix 9

Q. How to compute determinant ?

𝑎11 𝑎12 𝑎13
For example, determinant of 3×3 matrix: 𝐴 = 𝑎21 𝑎22 𝑎23 can be found as
𝑎31 𝑎32 𝑎33
𝑎22 𝑎23 𝑎21 𝑎23 𝑎21 𝑎22
𝑑𝑒𝑡𝐴 = 𝑎11 𝑎 𝑎33 − 𝑎12 𝑎 𝑎33 + 𝑎13 𝑎 𝑎32
32 31 31
= 𝑎11 𝑎22 𝑎33 − 𝑎22 𝑎32 − 𝑎12 𝑎21 𝑎33 − 𝑎31 𝑎23 + 𝑎13 𝑎21 𝑎32 − 𝑎31 𝑎22
= 𝑎11 𝑎22 𝑎33 − 𝑎11 𝑎23 𝑎32 − 𝑎12 𝑎21 𝑎33 + 𝑎12 𝑎31 𝑎23 + 𝑎13 𝑎21 𝑎32 − 𝑎13 𝑎31 𝑎22 .
Now let’s focus of the first two terms specifically

𝑎11 𝑎12 𝑎13

With either cases, each term of the product occupies a unique row and
𝐴 = 𝑎21 𝑎22 𝑎23
column. This is true of the other terms in the above expression also.
𝑎31 𝑎32 𝑎33

MEC E 390 -A1

Cramer’s rule: example 10

Based on this observation, we find…

a. If a matrix is diagonal, upper triangular or lower triangular, determinant is the product of the
diagonal entries (i.e. 𝑎11 𝑎22 𝑎33, all other terms are nulled out (multiplication with zeros)).
b. If an entire row or column consist of the zeros ⇒ 𝑑𝑒𝑡𝐴 = 0.
c. If two rows (or columns) are interchanged the determinant change sign.

Determinant may be evaluated in matlab using either the built in function ‘det’ or a personally
designed code (see, for example, Lecture7_determinant.m).
Example 1. Find solutions of the following linear equations using Cramer’s Rule
2𝑥1 + 𝑥2 = 1
2𝑥1 + 3𝑥2 = 2
2 1 1 𝑏 𝑎12 1 1
Ans. 𝐴 = ,𝑏= ⇒ 𝑑𝑒𝑡𝐴 = 2 × 3 − 2 × 1 = 4, Also we find that A1 = 1 = and
2 3 2 𝑏2 𝑎22 2 3
2 1
𝐴2 = . Thus 𝑑𝑒𝑡𝐴1 = 1 × 3 − 2 × 1 = 1 and 𝑑𝑒𝑡𝐴2 = 2 × 2 − 2 × 1 = 2. Consequently, we obtain
2 2
𝑑𝑒𝑡𝐴1 1 𝑑𝑒𝑡𝐴2 2
𝑥1 = = = 0.25 and 𝑥2 = = = 0.5
𝑑𝑒𝑡𝐴 4 𝑑𝑒𝑡𝐴 4

MEC E 390 -A1

Example continued 11

Example 2. Find solutions of the following linear equations using Cramer’s Rule (show 𝑥1 ONLY)
2𝑥1 − 𝑥2 + 𝑥3 = 4
4𝑥1 + 3𝑥2 − 𝑥3 = 2
3𝑥1 + 2𝑥2 − 2𝑥3 = 15

2 −1 1
Ans. 𝑑𝑒𝑡𝐴 = 4 3 −1 = 2 3 × −2 − −1 × 2 − −1 4 × −2 − −1 × 3 + 1 4 × 2 − 3 × 3
3 2 −2
= 2 −4 + −5 + −1 = −14
4 −1 1
𝑑𝑒𝑡𝐴1
To find 𝑥1 we require to compute: 𝑥1 = , where 𝐴1 = 2 3 −1
𝑑𝑒𝑡𝐴
15 2 −2
Thus we find 𝑑𝑒𝑡𝐴1 = −46
46
𝑥1 = − = 3.2857
−14
And similarly for others
Try matlab code: Lecture7_cramer.m

MEC E 390 -A1

Invers of matrixes 12

Solving a system of linear equations requires

𝐴𝒙 = 𝒃 ⇒ 𝒙 = 𝐴−1 𝒃
Therefore it is essential to determine the inverse of a matrix.
 Calculation of 𝐴−1
𝑎𝑑𝑗 𝐴
From math 102, we learned that 𝐴−1 = where 𝑎𝑑𝑗 𝐴 denotes the adjoint of 𝐴 so that
det 𝐴
T
𝐴11 𝐴21 ⋯ 𝐴𝑛1 𝐴11 𝐴12 ⋯ 𝐴1𝑛
𝐴 𝐴22 ⋯ 𝐴2𝑛 𝐴 𝐴22 ⋯ 𝐴2𝑛
𝑎𝑑𝑗 𝐴 = 12 = 21
⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮
𝐴1𝑛 𝐴2𝑛 ⋯ 𝐴𝑛𝑛 𝐴𝑛1 𝐴𝑛2 ⋯ 𝐴𝑛𝑛
Here the cofactors 𝐴11 , 𝐴21 , … 𝐴𝑛𝑛 are defined by
𝑖+𝑘 𝑀
𝐴𝑖𝑘 = −1 𝑖𝑘

where 𝑀𝑖𝑘 is the minor. i.e. the determinant of the (n-1) by (n-1) submatrix of 𝐴 obtained by
deleting the 𝑖 𝑡ℎ row and the 𝑖 𝑡ℎ column.

MEC E 390 -A1

Example continued 13

Example 3. Find the inverse of

2 0 −1
𝐴 = 3 −5 1
0 4 2
Using the result in Example 2., the determinant of 𝐴 can be obtained as 𝑑𝑒𝑡𝐴 = −40. Also to obtain
the adjoint, we have to calculate 𝐴11 , 𝐴21 etc…
−5 1 3 1 3 −5
𝐴11 = = −14, 𝐴12 = = − 6, 𝐴13 = = 12 𝑒𝑡𝑐 …
4 2 0 2 0 4
−14 −6 12 T −14 −4 −5
∴ 𝑎𝑑𝑗 𝐴 = −4 4 −8 ⇒ 𝑎𝑑𝑗 𝑎 = −6 4 −5
−5 −5 −10 12 −8 −10

𝑎𝑑𝑗 𝐴 0.35 0.1 0.125

∴ 𝐴−1 = = 0.15 −0.1 0.125
det 𝐴
−0.3 0.2 0.25
Note: we can always check our calculation by verifying that
𝐴𝐴−1 = 𝐴−1 𝐴 = 𝐼 (Identity matrix)
Check also : Lecture7_adjoint.m &Lecture7_inverse.m

MEC E 390 -A1