Transcript Matrix Methods - Civil and Environmental Engineering | SIU

MATRIX METHODS
SYSTEMS OF LINEAR EQUATIONS
ENGR 351
Numerical Methods for Engineers
Southern Illinois University Carbondale
College of Engineering
Dr. L.R. Chevalier
Dr. B.A. DeVantier
Copyright© 2003 by Lizette R. Chevalier and Bruce A. DeVantier
Permission is granted to students at Southern Illinois University at Carbondale
to make one copy of this material for use in the class ENGR 351, Numerical
Methods for Engineers. No other permission is granted.
All other rights are reserved. No part of this publication may be reproduced,
stored in a retrieval system, or transmitted, in any form or by any means,
electronic, mechanical, photocopying, recording, or otherwise, without
the prior written permission of the copyright owner.
Systems of Linear Algebraic
Equations
Specific Study Objectives
• Understand the graphic interpretation of
ill-conditioned systems and how it
relates to the determinant
• Be familiar with terminology: forward
elimination, back substitution, pivot
equations and pivot coefficient
Specific Study Objectives
• Apply matrix inversion to evaluate stimulusresponse computations in engineering
• Understand why the Gauss-Seidel method is
particularly well-suited for large sparse
systems of equations
• Know how to assess diagonal dominance of a
system of equations and how it relates to
whether the system can be solved with the
Gauss-Seidel method
Specific Study Objectives
• Understand the rationale behind
relaxation and how to apply this
technique
How to represent a system of linear
equations as a matrix
[A]{x} = {c}
where {x} and {c} are both column vectors
How to represent a system of linear
equations as a matrix
0.3 x1  0.52 x2  x3  0.01
0.5 x1  x2  1.9 x3  0.67
0.1x1  0.3 x2  0.5 x3  0.44
 A{ X }  {C}
0.3 0.52 1   x1    0.01
  

0.5

1
1.9  x2    0.67 


 0.1 0.3 0.5  x3   0.44
Practical application
• Consider a problem in structural
engineering
• Find the forces and reactions associated
with a statically determinant truss
90
30
hinge: transmits both
vertical and horizontal
forces at the surface
60
roller: transmits
vertical forces
1000 kg
1
F1
H2 2
90
F3
60
30
3
F2
V2
FREE BODY DIAGRAM
V3
F  0
F  0
H
v
Node 1
F1,V
30
F1
F1,H
60
F3


F

0


F
cos
30

F
cos
60
 F1,H
 H
1
3


F

0


F
sin
30

F
sin
60
 F1,V
V
1
3
 F1 cos 30  F3 cos 60  0
 F1 sin 30  F3 sin 60  1000
Node 2
F1
30
F2
H2
V2

F

0

H

F

F
cos
30
 H
2
2
1

F

0

V

F
sin
30
V
2
1
Node 3
F3
60
F2
V3

F

0


F
cos
60
 F2
 H
3

F

0

F
sin
60
 V3
V
3
 F1 cos 30  F3 cos 60  0
 F1 sin 30  F3 sin 60  1000
H 2  F2  F1 cos 30  0
V2  F1 sin 30  0
 F3 cos 60  F2  0
F3 sin 60  V3  0
SIX EQUATIONS
SIX UNKNOWNS
Do some book keeping
F1
F2
F3
H2
V2
V3
1
-cos30
0
cos60
0
0
0
0
2
-sin30
0
-sin60
0
0
0
-1000
3
cos30
1
0
1
0
0
0
4
sin30
0
0
0
1
0
0
5
0
-1
-cos60
0
0
0
0
6
0
0
sin60
0
0
1
0
This is the basis for your matrices and the equation
[A]{x}={c}
0.866
 0.5

 0.866
 0.5

 0
 0

0
0
0.5
0.866
0
0
0
0
0
0
1
0
0
0
1
0
0
1
0
0
1
0
0.5
0.866
0
0
0
0
0
1
  F1   0 
  F  1000
 2  

  F3   0 
 H    0 

 2  
 V2   0 
V   0 
 3  

System of Linear Equations
• We have focused our last lectures on
finding a value of x that satisfied a
single equation
• f(x) = 0
• Now we will deal with the case of
determining the values of x1, x2, .....xn,
that simultaneously satisfy a set of
equations
System of Linear Equations
• Simultaneous equations
•
•
•
•
f1(x1, x2, .....xn) = 0
f2(x1, x2, .....xn) = 0
.............
fn(x1, x2, .....xn) = 0
• Methods will be for linear equations
• a11x1 + a12x2 +...... a1nxn =c1
• a21x1 + a22x2 +...... a2nxn =c2
• ..........
• an1x1 + an2x2 +...... annxn =cn
Mathematical Background
Matrix Notation
 a11 a12
a
a22
21
 A  
.
 .
a
 m1 am2
•
•
•
•
a1n 
a23 ... a2 n 

.
. 
am3 ... amn 
a13 ...
a horizontal set of elements is called a row
a vertical set is called a column
first subscript refers to the row number
second subscript refers to column number
 a11 a12 a13 ... a1n 
a a

a
...
a
21
22
23
2n 

A 
 .
.
.
. 


a m1 a m 2 a m3 ... a mn 
This matrix has m rows an n column.
It has the dimensions m by n (m x n)
 a11 a12 a13 ... a1n 
 a a a ... a 
A    21 22 23 2n 
 .
.
.
. 


a m1 a m 2 a m 3 ... a mn 
This matrix has m rows and n column.
It has the dimensions m by n (m x n)
note
subscript
column 3
 a11 a12
a
a22
21
 A  
.
 .
a
 m1 am2
a1n 
a23 ... a2 n 

.
. 
am3 ... amn 
Note the consistent
scheme with subscripts
denoting row,column
a13 ...
row 2
Types of Matrices
Row vector: m=1
 B  b1
b2
....... bn 
Column vector: n=1
 c1 
c 
 2
C   . 
.
 
cm 
Square matrix: m = n
a11 a12
 A  a21 a22
a31 a32
a13 
a23 

a33 
Definitions
•
•
•
•
•
•
•
•
Symmetric matrix
Diagonal matrix
Identity matrix
Inverse of a matrix
Transpose of a matrix
Upper triangular matrix
Lower triangular matrix
Banded matrix
Symmetric Matrix
aij = aji for all i’s and j’s
5 1 2 
 A  1 3 7
2 7 8
Does a23 = a32 ?
Yes. Check the other elements
on your own.
Diagonal Matrix
A square matrix where all elements
off the main diagonal are zero
 a11
0
 A  
0
0

0
0
a22
0
0
0
a33
0
0
0

0
a44 
Identity Matrix
A diagonal matrix where all elements
on the main diagonal are equal to 1
1
0
 A  
0
0

0 0 0
1 0 0

0 1 0
0 0 1
The symbol [I] is used to denote the identify matrix.
Inverse of [A]
 A A
1
  A  A  I 
1
Transpose of [A]
 a11
a
 12
 .
t
 A  
.

 .

 a1n
a21
.
.
.
a22
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
a2 n
.
.
.
am1 
am2 

. 
. 

. 

amn 
Upper Triangle Matrix
Elements below the main diagonal
are zero
a11 a12
 A   0 a22
0
 0
a13 
a23 

a33 
Lower Triangular Matrix
All elements above the main diagonal
are zero
5 0 0
 A  1 3 0
2 7 8
Banded Matrix
All elements are zero with the
exception of a band centered on the
main diagonal
 a11 a12
a
a22
21
 A  
 0 a32
0
0

0
a23
a33
a43
0
0

a34 
a44 
Matrix Operating Rules
• Addition/subtraction
• add/subtract corresponding terms
• aij + bij = cij
• Addition/subtraction are commutative
• [A] + [B] = [B] + [A]
• Addition/subtraction are associative
• [A] + ([B]+[C]) = ([A] +[B]) + [C]
Matrix Operating Rules
• Multiplication of a matrix [A] by a scalar
g is obtained by multiplying every
element of [A] by g
 ga11
 ga
 21
 .
 B  g A  
.

 .

 gam1
ga12
.
.
.
ga22
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
gam2
.
.
.
ga 1n 
ga2 n 

. 
. 

. 

gamn 
Matrix Operating Rules
• The product of two matrices is represented as
[C] = [A][B]
• n = column dimensions of [A]
• n = row dimensions of [B]
N
c ij   a ik b kj
k 1
Simple way to check whether
matrix multiplication is possible
exterior dimensions conform to dimension of resulting matrix
[A] m x n [B] n x k = [C] m x k
interior dimensions
must be equal
Recall the equation presented for
matrix multiplication
• The product of two matrices is represented as
[C] = [A][B]
• n = column dimensions of [A]
• n = row dimensions of [B]
N
c ij   a ik b kj
k 1
Example
Determine [C] given [A][B] = [C]
 1

A    1
 0
 2

B   3
 3
2

4 2
 2 3
4 1

2  1
0 2 
3
Matrix multiplication
• If the dimensions are suitable, matrix
multiplication is associative
• ([A][B])[C] = [A]([B][C])
• If the dimensions are suitable, matrix
multiplication is distributive
• ([A] + [B])[C] = [A][C] + [B][C]
• Multiplication is generally not
commutative
• [A][B] is not equal to [B][A]
Determinants
Denoted as det A or lAl
for a 2 x 2 matrix
a b
c d
a b
c d
 ad  bc
 ad  bc
Determinants
For a 3 x 3
+
2
3
1
9 2
3
4
5 2
2
+
3
1
9 2
3
4
5 2
2
3
1
9 2
3
4
5 2
Problem
1 7 9 
 4 3 2


6 1 5
Determine the determinant of the matrix.
Properties of Determinants
• det A = det AT
• If all entries of any row or column is
zero, then det A = 0
• If two rows or two columns are
identical, then det A = 0
• Note: determinants can be calculated
using mdeterm function in Excel
Excel Demonstration
• Excel treats matrices as arrays
• To obtain the results of multiplication,
addition, and inverse operations, you hit
control-shift-enter as opposed to enter.
• The resulting matrix cannot be
altered…let’s see an example using
Excel in class
matrix.xls
Matrix Methods
•
•
•
•
Cramer’s Rule
Gauss elimination
Matrix inversion
Gauss Seidel/Jacobi
a11 a12
 A  a21 a22
a31 a32
a13 
a23 

a33 
Graphical Method
2 equations, 2 unknowns
a11 x1  a12 x2  c1
a21 x1  a22 x2  c2
 a11 
c1
 x1 
x2   
a12
 a12 
 a21 
c2
 x1 
x2   
a22
 a22 
x2
( x1, x2 )
x1
x2
3x1  2 x2  18
 3
x2    x1  9
 2
9
3
2
x1
x2
 x1  2 x2  2
 1
x2    x1  1
 2
1
2
1
x1
x2
3x1  2 x2  18
 x1  2 x2  2
9
3
 3
x2    x1  9
 2
 1
x2    x1  1
 2
2
(4,3)
1
2
1
x1
Check: 3(4) + 2(3) = 12 + 6 = 18
Special Cases
• No solution
• Infinite solution
• Ill-conditioned
x2
( x1, x2 )
x1
a) No solution - same slope
b) infinite solution
f(x)
x
f(x)
-1/2 x1 + x2 = 1
-x1 +2x2 = 2
c) ill conditioned
so close that the points of
intersection are difficult to
detect visually
f(x)
x
x
Ill Conditioned Systems
• If the determinant is zero, the slopes
are identical
a 11 x 1  a 12 x 2  c1
a 21 x 1  a 22 x 2  c2
Rearrange these equations so that we have an
alternative version in the form of a straight line:
i.e. x2 = (slope) x1 + intercept
Ill Conditioned Systems
a11
c1
x2  
x1 
a12
a12
x2  
a21
c
x1  2
a22
a22
If the slopes are nearly equal (ill-conditioned)
a11 a21

a12 a22
a11a22  a21a12
a11a22  a21a12  0
Isn’t this the determinant?
a11
a12
a21
a22
 det A
Ill Conditioned Systems
If the determinant is zero the slopes are equal.
This can mean:
- no solution
- infinite number of solutions
If the determinant is close to zero, the system is ill
conditioned.
So it seems that we should use check the determinant of a
system before any further calculations are done.
Let’s try an example.
Example
Determine whether the following matrix is ill-conditioned.
37.2 4.7  x1   22 
19.2 2.5  x    12

 2  

Solution
37.2 4.7
19.2
2.5
 37.22.5  4.719.2
 2.76
What does this tell us? Is this close to zero? Hard to say.
If we scale the matrix first, i.e. divide by the largest
a value in each row, we can get a better sense of things.
Solution
x
0
5
10
15
0
y
-20
-40
-60
This is further justified
when we consider a graph
of the two functions.
-80
Clearly the slopes are
nearly equal
1 0126
.
 0.004
1 0130
.
Another Check
• Scale the matrix of coefficients, [A], so that the
largest element in each row is 1. If there are
elements of [A]-1 that are several orders of
magnitude greater than one, it is likely that the
system is ill-conditioned.
• Multiply the inverse by the original coefficient matrix.
If the results are not close to the identity matrix, the
system is ill-conditioned.
• Invert the inverted matrix. If it is not close to the
original coefficient matrix, the system is illconditioned.
We will consider how to obtain an inverted matrix later.
Cramer’s Rule
• Not efficient for solving large numbers
of linear equations
• Useful for explaining some inherent
problems associated with solving linear
equations.
 a11
a
 21
a31
a12
a22
a32
a13   x1  b1 
   

a23  x2   b2   Ax  b

  
a33  
 x3  b3 
Cramer’s Rule
b1
1
x1 
b2
A
b3
a12
a 22
a13
a 23
a 32
a 33
to solve for
xi - place {b} in
the ith column
Cramer’s Rule
b1
x1 
1
b2
A
b3
a11
1
x3  a21
A
a31
a12
a22
a32
a12
a22
a32
a13
a11
b1
1
a23 x2  a21 b2
A
a33
a31 b3
b1
b2
b3
a13
a23
a33
to solve for
xi - place {b} in
the ith column
EXAMPLE
Use of Cramer’s Rule
2 x1  3x2  5
x1  x2  5
2 3  x1  5
1 1   x   5

 2   
Solution
2 3  x1  5
1 1   x   5

 2   
A  21   31  2  3  5
1 5 3  1
20
x1 
   51   35 
4
5 5 1  5
5
1 2 5  1
5
x2 
   25  51   1
5 1 5  5
5
Elimination of Unknowns
( algebraic approach)
a11 x1  a12 x2  c1
a21 x1  a22 x2  c2
a11 x1  a12 x2  c1 a21 
a21 x1  a22 x2  c2 a11 
a21a11 x1  a21a12 x2  a21c1 SUBTRACT
a21a11 x1  a11a22 x2  a11c2
Elimination of Unknowns
( algebraic approach)
a21a11 x1  a21a12 x2  a21c1 SUBTRACT
a21a11 x1  a11a22 x2  a11c2
a11a21 x2  a22 a11 x2  c1a21  c2 a11
a21c1  a11c2
x2 
a12 a21  a22 a11
x1 
a22 c1  a12 c2
a11a22  a12 a21
NOTE: same result as
Cramer’s Rule
Gauss Elimination
• One of the earliest methods developed
for solving simultaneous equations
• Important algorithm in use today
• Involves combining equations in order
to eliminate unknowns and create an
upper triangular matrix
• Progressively back substitute to find
each unknown
Two Phases of Gauss
Elimination
 a11
a
 21
 a31
 a11
0

 0
a12
a22
a32
a12
'
a22
0
a13 | c1 
a23 | c2 

a33 | c3 
a13 | c1 
'
| c2' 
a23

''
''
a33 | c3 
Forward
Elimination
Note: the prime indicates
the number of times the
element has changed from
the original value.
Two Phases of Gauss
Elimination
a11
0

 0
a12
'
a22
0
c3''
x3  ''
a33
x2

c

x1 
'
2
a13 | c1 
'
a23
| c2' 

''
''
a33 | c3 
 a123 x3 
'
a22
c1  a12 x2  a13 x3 
a11
Back substitution
Rules
• Any equation can be multiplied (or
divided) by a nonzero scalar
• Any equation can be added to (or
subtracted from) another equation
• The positions of any two equations in
the set can be interchanged.
EXAMPLE
Perform Gauss Elimination of the following matrix.
2 x1  x2  3x3  1
4 x1  4 x2  7 x3  1
2 x1  5x2  9 x3  3
Solution
2 x1  x2  3x3  1
4 x1  4 x2  7 x3  1
2 x1  5x2  9 x3  3
Multiply the first equation by
a21/ a11 = 4/2 = 2
Note: a11 is called the pivot element
4x1  2x2  6x3  2
a21 / a11 = 4/2 = 2
2 x1  x2  3x3  1
4 x1  4 x2  7 x3  1
2 x1  5x2  9 x3  3
4 x1  2 x2  6 x3  2
a21 / a11 = 4/2 = 2
2 x1  x2  3x3  1
4 x1  2 x 2  6 x3  2
4 x1  4 x2  7 x3  1
4 x1  4 x 2  7 x3  1
2 x1  5x2  9 x3  3
2 x1  5 x 2  9 x3  3
a21 / a11 = 4/2 = 2
2 x1  x2  3x3  1
4 x1  2 x 2  6 x3  2
4 x1  4 x2  7 x3  1
4 x1  4 x 2  7 x3  1
2 x1  5x2  9 x3  3
2 x1  5 x 2  9 x3  3
Subtract the revised first equation from the
second equation
a21 / a11 = 4/2 = 2
2 x1  x2  3x3  1
4 x1  2 x 2  6 x3  2
4 x1  4 x2  7 x3  1
4 x1  4 x 2  7 x3  1
2 x1  5x2  9 x3  3
2 x1  5 x 2  9 x3  3
Subtract the revised first equation from the
second equation
4  4 x1  4  2 x2  7  6 x3  1  2
0x1  2 x2  x3  1
a21 / a11 = 4/2 = 2
2 x1  x2  3x3  1
4 x1  2 x 2  6 x3  2
4 x1  4 x2  7 x3  1
4 x1  4 x 2  7 x3  1
2 x1  5x2  9 x3  3
2 x1  5 x 2  9 x3  3
Subtract the revised first equation from the
second equation
4  4 x1  4  2 x2  7  6 x3  1  2
0x1  2 x2  x3  1
NEW
MATRIX
2 x1  x2  3x3  1
0 x1  2 x2  x3  1
2 x1  5 x2  9 x3  3
a21 / a11 = 4/2 = 2
2 x1  x2  3x3  1
4 x1  2 x 2  6 x3  2
4 x1  4 x2  7 x3  1
4 x1  4 x 2  7 x3  1
2 x1  5x2  9 x3  3
2 x1  5 x 2  9 x3  3
Subtract the revised first equation from the
second equation
4  4 x1  4  2 x2  7  6 x3  1  2
0x1  2 x2  x3  1
NOW LET’S
GET A ZERO
HERE
2 x1  x2  3x3  1
0 x1  2 x2  x3  1
2 x1  5 x 2  9 x3  3
Solution
2 x1  x2  3x3  1
4 x1  4 x2  7 x3  1
2 x1  5x2  9 x3  3
Multiply equation 1 by a31/a11 = 2/2 = 1
and subtract from equation 3
2  2 x1  5  1 x2  9  3 x3  3  1
0x1  4 x2  6 x3  2
Solution
Following the same rationale, subtract the 3rd equation
from the first equation
2 x1  x2  3x3  1
4 x1  4 x2  7 x3  1
2 x1  5x2  9 x3  3
2 x1  x2  3x3  1
2 x 2  x 3  1
4 x2  6x3  2
Continue the
computation
by multiplying
the second equation
by a32’/a22’ = 4/2 =2
Subtract the third
equation of the new
matrix
Solution
2 x1  x2  3x3  1
2 x 2  x 3  1
4 x2  6x3  2
2 x1  x2  3x3  1
2 x 2  x 3  1
4 x3  4
THIS DERIVATION OF
AN UPPER TRIANGULAR MATRIX
IS CALLED THE FORWARD
ELIMINATION PROCESS
Solution
From the system we immediately calculate:
2 x1  x2  3x3  1
2 x2  x 3  1
4 x3  4
4
x3   1
4
Continue to back substitute
1  1
x2 
 1
2
1  3   1
1
x1 

2
2
THIS SERIES OF
STEPS IS THE
BACK
SUBSTITUTION
Pitfalls of the Elimination Method
• Division by zero
• Round off errors
• magnitude of the pivot element is small compared
to other elements
• Ill conditioned systems
Pivoting
• Partial pivoting
• rows are switched so that the pivot element is not
zero
• rows are switched so that the largest element is
the pivot element
• Complete pivoting
• columns as well as rows are searched for the
largest element and switched
• rarely used because switching columns changes
the order of the x’s adding unjustified complexity
to the computer program
For example
Pivoting is used here to avoid division
by zero
2 x2  3x3  8
4 x1  6 x2  7 x3  3
2 x1 
x2  6 x3  5
4 x1  6 x2  7 x3  3
2 x2  3 x3  8
2 x1  x2  6 x3  5
Another Improvement: Scaling
• Minimizes round-off errors for cases where
some of the equations in a system have much
larger coefficients than others
• In engineering practice, this is often due to
the widely different units used in the
development of the simultaneous equations
• As long as each equation is consistent, the
system will be technically correct and solvable
Example
(solution in notes)
Use Gauss Elimination to solve the following set
of linear equations. Employ partial pivoting when
necessary.
3x2  13x3  50
2 x1  6x2  x3  45
4 x1  8x3  4
Solution
3x2  13x3  50
2 x1  6x2  x3  45
4 x1  8x3  4
First write in matrix form, employing short hand
presented in class.
 0
 2

 4
3
13

6
1

0
8

50
45 

4 
We will clearly run into
problems of division
by zero.
Use partial pivoting
 0
 2

 4
3
6
13
1


0
8

50
45 

4 
Pivot with equation
with largest an1
 0
 2

 4
 4
 2

 0
3
 13

6
1

0
8

0
8

6
1

3
 13

 50
45 
4 
4 
45 
 50
 0
 2

 4
 4
 2

 0
 4
 0

 0
3
 13

6
0
1
8


0
8

6
3
1
 13


0
6
8
3


3
 13

 50
45 
4 
4 
45 
 50
4 
43 
 50
Begin developing
upper triangular matrix
 4
 0

 0
 4
 0

 0
0
8

6
3
3
13


4 
43 

50
0
8

6
0
3
14 .5


4 
43 

28.5
43  31.966
28.5
 1.966 x2 
 8149
.
14 .5
6
4  81.966
x1 
 2 .931
4
CHECK
x3 
3 8149
.   131.966  50 okay
...end of
problem
Gauss-Jordan
• Variation of Gauss elimination
• Primary motive for introducing this method is
that it provides a simple and convenient
method for computing the matrix inverse.
• When an unknown is eliminated, it is
eliminated from all other equations, rather
than just the subsequent one
Gauss-Jordan
a11 a12
 A   0 a22
0
 0
a13 
a23 

a33 
1
0
 A  
0
0

0 0 0
1 0 0

0 1 0
0 0 1
• All rows are normalized by dividing them by
their pivot elements
• Elimination step results in an identity matrix
rather than an UT matrix
Graphical depiction of Gauss-Jordan
a11 a12
'
a
a
 21 22
a31 a32
a13 | c1 
'
a23
| c2' 

''
''
a33 | c3 
1 0 0 | c1n  

 n 
0
1
0
|
c
2 

0 0 1 | c3 n  
Graphical depiction of Gauss-Jordan
a11 a12
'
a
a
 21 22
a31 a32
a13 | c1 
'
a23
| c2' 

''
''
a33 | c3 
1 0 0 | c1n  

 n 
0
1
0
|
c
2 

0 0 1 | c3 n  
1 0 0 | c1n  

n 
0
1
0
|
c
2 

0 0 1 | c3 n  
 c1n 
x1
x2
 c2 n 
x3  c3 n 
Matrix Inversion
• [A] [A] -1 = [A]-1 [A] = I
• One application of the inverse is to solve
several systems differing only by {c}
• [A]{x} = {c}
• [A]-1[A] {x} = [A]-1{c}
• [I]{x}={x}= [A]-1{c}
• One quick method to compute the inverse is
to augment [A] with [I] instead of {c}
Graphical Depiction of the Gauss-Jordan
Method with Matrix Inversion
 A
 a11
a
 21
a31
 1

 0
 0
I 
a12
a13

1
0
a22
a23

0
1
a32
a33

0
0
0
0

a111
a121
1
0

1
a21
1
a22
0
1

1
1
a31
a 32
I 
 A1
0
0

1 
Note: the superscript
“-1” denotes that
the original values
a131  have been converted
to the matrix inverse,
1 
a23  not 1/a
ij
1
a33 
WHEN IS THE INVERSE MATRIX USEFUL?
CONSIDER STIMULUS-RESPONSE
CALCULATIONS THAT ARE SO COMMON IN
ENGINEERING.
Stimulus-Response Computations
• Conservation Laws
mass
force
heat
momentum
• We considered the conservation
of force in the earlier example of
a truss
Stimulus-Response Computations
• [A]{x}={c}
• [interactions]{response}={stimuli}
• Superposition
• if a system subject to several different stimuli, the
response can be computed individually and the
results summed to obtain a total response
• Proportionality
• multiplying the stimuli by a quantity results in the
response to those stimuli being multiplied by the
same quantity
• These concepts are inherent in the scaling of terms
during the inversion of the matrix
Example
Given the following, determine {x} for the
two different loads {c}
Ax  c
 2 1 1 
A 1   2 6
3 
  3 1  4
c  1
T
c  4
T
2 3
 7 1
Solution
Ax  c
 2 1 1 


1
A   2 6
3
  3 1  4
c  1
cT  4
T
2 3
 7 1
{c}T = {1 2 3}
x1 = (2)(1) + (-1)(2) + (1)(3) = 3
x2 = (-2)(1) + (6)(2) + (3)(3) = 19
x3 = (-3)(1) + (1)(2) + (-4)(3) = -13
{c} T= {4 -7 1)
x1 = (2)(4) + (-1)(-7) + (1)(1)=16
x2 = (-2)(4) + (6)(-7) + (3)(1) = -47
x3 = (-3)(4) + (1)(-7) + (-4)(1) = -23
Gauss Seidel Method
• An iterative approach
• Continue until we converge within some prespecified tolerance of error
• Round off is no longer an issue, since you control
the level of error that is acceptable
• Fundamentally different from Gauss elimination.
This is an approximate, iterative method
particularly good for large number of equations
Gauss-Seidel Method
• If the diagonal elements are all nonzero, the
first equation can be solved for x1
c1  a12 x2  a13x3  a1n xn
x1 
a11
• Solve the second equation for x2, etc.
To assure that you understand this, write the equation for x2
c1  a12 x2  a13 x 3  a1n x n
x1 
a11
c2  a 21 x1  a 23 x3  a 2 n xn
x2 
a 22
c3  a 31 x1  a 32 x 2  a 3n xn
x3 
a 33

cn  a n1 x1  a n 3 x2  a nn 1 x n 1
xn 
a nn
Gauss-Seidel Method
• Start the solution process by guessing
values of x
• A simple way to obtain initial guesses is
to assume that they are all zero
• Calculate new values of xi starting with
• x1 = c1/a11
• Progressively substitute through the
equations
• Repeat until tolerance is reached
Gauss-Seidel Method
x1  c1  a12 x2  a13 x3  / a11
x2  c2  a21 x1  a23 x3  / a22
x3  c3  a31 x1  a32 x2  / a33
x1  c1  a12 0  a13 0 / a11 
c1
a11
x2  c2  a21 x '1  a23 0 / a22  x '2
x3  c3  a31 x '1  a32 x '2  / a33  x '3
 x '1
Example
Given the following augmented matrix,
complete one iteration of the Gauss Seidel
method.
2
4

 3
3
1
2
1
2
1



2
2 

1 
2
4

 3
3
1
2
1
2
1



2
2

1 
GAUSS SEIDEL
x1  c1  a12 0  a13 0 / a11 
c1
a11
 x '1
2   30  10 2
x1 
 1
2
2
x2  c2  a21 x '1  a23 0 / a22  x '2
2   41   20 2  4
x2 

 6
1
1
x3  c3  a31 x '1  a32 x '2  / a33  x '3
1   31   2 6 1  3  12
x3 

 10
1
1
Jacobi Iteration
• Iterative like Gauss Seidel
• Gauss-Seidel immediately uses the
value of xi in the next equation to
predict x i+1
• Jacobi calculates all new values of xi’s to
calculate a set of new xi values
Graphical depiction of difference between Gauss-Seidel and Jacobi
FIRST ITERATION
x1   c1  a12 x2  a13x3  / a11
x1   c1  a12 x2  a13x3  / a11
x2   c2  a 21x1  a 23x3  / a 22
x2   c2  a 21x1  a 23x3  / a 22
x3   c3  a 31x1  a 32 x2  / a 33
x3   c3  a 31x1  a 32 x2  / a 33
SECOND ITERATION
x1   c1  a12 x2  a13x3  / a11
x1   c1  a12 x2  a13x3  / a11
x2   c2  a 21x1  a 23x3  / a 22
x2   c2  a 21x1  a 23x3  / a 22
x3   c3  a 31x1  a 32 x2  / a 33
x3   c3  a 31x1  a 32 x2  / a 33
Example
Given the following augmented matrix, complete
one iteration of the Gauss Seidel method and
the Jacobi method.
2
4

 3
3
1
2
1
2
1



2
2 

1 
Note: We worked the Gauss Seidel method earlier
Gauss-Seidel Method
convergence criterion
 a ,i
xij  xij 1

 100   s
j
xi
as in previous iterative procedures in finding the roots,
we consider the present and previous estimates.
As with the open methods we studied previously with one
point iterations
1. The method can diverge
2. May converge very slowly
Convergence criteria for two
linear equations
c1 a12
u x1 , x2  

x2
a11 a11
c2 a21
v x1 , x2  

x2
a22 a22
Class question:
where do these
formulas come from?
consider the partial derivatives of u and v
u
u
a12
0

x1
x 2
a11
v
a21

x1
a22
v
0
x 2
Convergence criteria for two linear
equations cont.
u v

1
x x
u v

1
y y
Criteria for convergence
where presented earlier
in class material
for nonlinear equations.
Noting that x = x1 and
y = x2
Substituting the previous equation:
Convergence criteria for two linear
equations cont.
a21
1
a22
a12
1
a11
This is stating that the absolute values of the slopes must
be less than unity to ensure convergence.
Extended to n equations:
aii   aij where j  1, n excluding j  i
Convergence criteria for two linear
equations cont.
aii   aij where j  1, n excluding j  i
This condition is sufficient but not necessary; for convergence.
When met, the matrix is said to be diagonally dominant.
Diagonal Dominance
  1 0.2 0.4   x1   9
 0. 2 0. 8 0. 4   x    3 

 2   
 0.1 0.5  0.9  x 3   4 
To determine whether a matrix is diagonally
dominant you need to evaluate the values on
the diagonal.
Diagonal Dominance
  1 0.2 0.4   x1   9
 0. 2 0. 8 0. 4   x    3 

 2   
 0.1 0.5  0.9  x 3   4 
Now, check to see if these numbers satisfy the
following rule for each row (note: each row
represents a unique equation).
aii   aij where j  1, n excluding j  i
Review the concepts
of divergence and
convergence by graphically
illustrating Gauss-Seidel
for two linear equations
x2
x1
u: 11x1  13x2  286
v: 11x1  9 x2  99
v: 11x1  9 x2  99
u: 11x1  13x2  286
x2
Note: we are converging
on the solution
x1
CONVERGENCE
u: 11x1  13x2  286
x2
v: 11x1  9 x2  99
Change the order of
the equations: i.e. change
direction of initial
estimates
x1
DIVERGENCE
Improvement of Convergence
Using Relaxation
This is a modification that will enhance slow convergence.
After each new value of x is computed, calculate a new value
based on a weighted average of the present and previous
iteration.
xinew  xinew  1    xiold
Improvement of Convergence Using
Relaxation
xinew  xinew  1    xiold
• if  = 1unmodified
• if 0 <  < 1
underrelaxation
• nonconvergent systems may converge
• hasten convergence by dampening out
oscillations
• if 1<  < 2
overrelaxation
• extra weight is placed on the present value
• assumption that new value is moving to the
correct solution by too slowly