Linear Equations
PSCI 702
September 21, 2005
Numerical Errors
• Round off Error
• Truncation Error
Linear Algebraic Equations
a11 x1 a12 x2 a13 x3 ... a1n xn c1
a21 x1 a22 x1 a23 x3 ... a2 n xn c2
a31 x1 a32 x1 a33 x3 ... a3n xn c3
an1 x1 an 2 x2 an 3 x3 ... ann xn cn
a11 a12 a13 a1n x1 c1
a a22 a23 a2 n x2 c 2
21
a31 a32 a33 a3n x3 c 3
an1
an 2 an 3 ann xn cn
Linear Algebraic Equations
• In Matrix Format
Ax c
• Solution:
1 1 1
A Ax A c x A c
Cramer’s Rule
• Mij is the determinant of the Matrix A with
the ith row and jth column removed.
• (-1)ij is called the cofactor of element aij.
n
Det A (1) aij M ij , j
i j
i 1
n n
Det A C ij aij , j aij C ij , i
i 1 i j
Cramer’s Rule
a11 a12 a13 a1n a11x1 a12 a13 a1n
a21 a22 a23 a2 n a21x1 a22 a23 a2 n
x1 a31 a32 a33 a3n a31x1 a32 a33 a3n
an1 an 2 an3 ann an1x1 an 2 an3 ann
a11 x1 a12 x2 a13 x3 ... a1n xn a12 a1n c1 a12 a13 a1n
a21 x1 a22 x1 a23 x3 ... a2 n xn a22 a2 n c2 a22 a23 a2 n
a31 x1 a32 x1 a33 x3 ... a3n xn a32 a3n c3 a32 a33 a3n
an1 x1 an 2 x2 an3 x3 ... ann xn an 2 ann cn an 2 an3 ann
Cramer’s Rule
c1 a12 a13 a1n a11 a1 j 1 c1 a1 j 1 a1n
c2 a22 a23 a2 n a21 a2 j 1 c2 a2 j 1 a2 n
c3 a32 a33 a3n a31 a3 j 1 c3 a3 j 1 a3n
cn an 2 an 3 ann an1 anj 1 cn an j 1 ann
x1 ,xj
a11 a12 a13 a1n a11 a12 a13 a1n
a21 a22 a23 a2 n a21 a22 a23 a2 n
a31 a32 a33 a3n a31 a32 a33 a3n
an1 an 2 an3 ann an1 an 2 an 3 ann
Cramer’s Rule
• 3n2 operations for the determinant.
• 3n3 operations for every unknown.
• Unstable for large Matrices.
• Large error propagation.
• Good for small Matrices (n<20).
Gaussian Elimination
• Divide each row by the leading element.
• Subtract row 1 from all other rows.
• Move to the second row an continue the
process.
a11 a12 a13 a1n c1
a a22 a23 a2 n c 2
21
a31 a32 a33 a3n c 3
an1 an 2 an 3 ann cn
Gaussian Elimination
a11 a12 a13 a11 c1
a a
a11 a11 a11
11 11
a31 a32 a33 a3n c3
a31 a31 a31 a31 a31
a an 2 an 3 ann cn
n1
an1
an1 an1 an1 an1
Gaussian Elimination
a12 a13 a11 c1
1 a
a11 a11 a11
11
0 a32 a12 a33 a13 a3n a11 c3 c1
a31 a11 a31 a11 a31 a11 a31 a11
an 2 a12 an 3 a13 ann a11 cn c1
0
an1 a11 an1 a11 an1 a11 an1 a11
Gaussian Elimination
1 a12 a13 a1n c1
0 1
a23 a2 n c2
0 0 1 a3n c3
0 0
0 1 cn
• Back substitution
x n cn
n
x i ci a x
j i 1
ij j
Gaussian Elimination
• Division by zero: May occur in the forward
elimination steps.
• Round-off error: Prone to round-off errors.
Gaussian Elimination
Consider the system of equations:
Use five significant figures with chopping
10 7 0 x1 7
3 2.099 6
x2 = 3.901
5
1 5 x3 6
At the end of Forward Elimination
10 7 0 x1 7
0 0.001 6 x = 6.001
2
0
0 15005
x3 15004
Gaussian Elimination
Back Substitution
15004
x3 0.99993
7 15005
10 0 x1 7
0 0.001 6 x 2 6.001
6.001 6 x3
0 15005 x3 15004 x2 1.5
0 0.001
7 7 x 2 0 x3
x1 0.3500
10
Gaussian Elimination
Compare the calculated values with the exact solution
x1 0
X exact x 1
2
x3 1
x1 0.35
X calculated x 2 1.5
x3 0.99993
Improvements
Increase the number of significant digits
Decreases round off error
Does not avoid division by zero
Gaussian Elimination with Partial Pivoting
Avoids division by zero
Reduces round off error
Partial Pivoting
Gaussian Elimination with partial pivoting applies row switching to
normal Gaussian Elimination.
How?
At the beginning of the kth step of forward elimination, find the maximum of
.......,ank
akk , ak 1,k ,.........
If the maximum of the values is a pk In the pth row, k p n,
then switch rows p and k.
Partial Pivoting
What does it Mean?
Gaussian Elimination with Partial Pivoting ensures that
each step of Forward Elimination is performed with the
pivoting element |akk| having the largest absolute value.
Partial Pivoting: Example
Consider the system of equations
10x1 7 x 2 7
3x1 2.099x 2 3x 3 3.901
5x 1 x 2 5x 3 6
In matrix form
10 7 0 x1 7
3 2.099 6 x 3.901
2 =
5
1 5 x3
6
Solve using Gaussian Elimination with Partial Pivoting using five
significant digits with chopping
Partial Pivoting: Example
Forward Elimination: Step 1
Examining the values of the first column
|10|, |-3|, and |5| or 10, 3, and 5
The largest absolute value is 10, which means, to follow the
rules of Partial Pivoting, we switch row1 with row1.
Performing Forward Elimination
10 7 0 x1 7 10 7 0 x1 7
3 2.099 6 x 3.901
5
2
1 5 x3 6
0 0.001 6 x 6.001
0
2.5
2
5 x3 2.5
Partial Pivoting: Example
Forward Elimination: Step 2
Examining the values of the first column
|-0.001| and |2.5| or 0.0001 and 2.5
The largest absolute value is 2.5, so row 2 is switched with
row 3
Performing the row swap
10 7 0 x1 7 10 7 0 x1 7
0 0.001 6 x 6.001
0
2
5 x3 2.5
0
2.5 5 x2 2.5
0 0.001 6 x3 6.001
2.5
Partial Pivoting: Example
Forward Elimination: Step 2
Performing the Forward Elimination results in:
10 7 0 x1 7
0 2.5 5 x 2 2.5
0
0 6.002 x3 6.002
Partial Pivoting: Example
Back Substitution
Solving the equations through back substitution
10 7
6.002
0 x1 7 x3 1
0 2.5 5 x 2 2.5 6.002
0
0 6.002 x3 6.002
2.5 5 x 2
x2 1
2.5
7 7 x 2 0 x3
x1 0
10
Partial Pivoting: Example
Compare the calculated and exact solution
The fact that they are equal is coincidence, but it does
illustrate the advantage of Partial Pivoting
x1 0 x1 0
X calculated x2 1
X exact x 2 1
x3 1
x3 1
Gauss Jordan Elimination
• Start with the following system of Matrices.
a11 a12 a13 a1n c1 1 0 0 0
a a22 a23 a2 n c 2 0 1 0 0
21
a31 a32 a33 a3n c 3 0 0 1 0
an1
an 2 an 3 ann cn 0
0 0 1
• Divide each row by the leading element.
• Subtract row 1 from all other rows.
• In addition to subtracting the line whose diagonal
term has been made unity from all those bellow it,
also subtract from the equations above it as well.
Gauss Jordan Elimination
a11 0 0 0 c1 b11 b12 b13 b1n
0 a 0 0 c2 b21 b22 b23 b2 n
22
0 0 a33 0 c3 b31 b32 b33 b3n
0
0 0 ann cn bn1
bn 2 bn 3 bnn
x1 c1
x c
2 2
x3 c3
xn cn
Matrix Factorization
• Assume A can be written as A=VU where V
and U are triangular Matrices.
a11 a12 a13 a1n v11 0 0 0 u11 u12 u13 u1n
a a2 n v21 v22 0 0 u22 u2 n
21 a22 a23 0 u23
a31 a32 a33 a3n v31 v32 v33 0 0 0 u33 u3 n
an1 an 2
an 3 ann vn1 vn 2
vn 3 vnn 0
0 0 unn
Matrix Factorization
Proof
If solving a set of linear equations AX C
If A V U Then V U X C
Multiply by V 1
Which gives V 1 V U X V 1 C
Remember V V I which leads to I U X V 1 C
1
Now, if I U U then U X V 1 C
Now, let V 1C Z
Which ends with V Z C (1)
and U X Z (2)
Matrix Factorization
How can this be used?
GivenAX C
Decompose A into V and U
Then solve VZ C for Z
And then solve U X Z for X
Matrix Factorization
How is this better or faster than Gauss
Elimination?
Let’s look at computational time.
n = number of equations
n3
To decompose [A], time is proportional to
3
To solve U X C and V Z C
n2
time is proportional to
2
Matrix Factorization
Therefore, total computational time for LU Decomposition is
proportional to 3 2
n3
n n
2( ) or n2
3 2 3
Gauss Elimination computation time is proportional to
3 2
n n
3 2
How is this better?
Matrix Factorization
What about a situation where the [C] vector changes?
In VU factorization, VU decomposition of [A] is independent of
the [C] vector, therefore it only needs to be done once.
Let m = the number of times the [C] vector changes
The computational times are proportional to
n3 n2 n3
VU factorization = m( ) Gauss Elimination= m( n 2 )
3 2 3
Consider a 100 equation set with 50 right hand side vectors
VU factorization = 8.33105 Gauss Elimination = 1.69 107
Matrix Factorization
Another Advantage
Finding the Inverse of a Matrix
VU Factorization Gauss Elimination
n3
4n 3 n3 n2 n4 n3
n( n )
2
n
3 3 2
3 3 2
For large values of n
4 3 3
n n 4n
3 2 3
Matrix Factorization
Method: Decompose [A] to [V] and [U]
1 0 0 u11 u12 u13
A V U v21 1
0 0 u22
u23
v31 v32
1 0
0 u33
[U] is the same as the coefficient matrix at the end of the forward
elimination step.
[V] is obtained using the multipliers that were used in the forward
elimination process
Matrix Factorization
Finding the [U] matrix
Using the Forward Elimination Procedure of Gauss Elimination
25 5 1
64 8 1
144 12 1
25 5 1
Row1
Row 2 (64) 0 4.8 1.56
25
144 12
1
25 5 1
Row1
Row 3 (144) 0 4.8 1.56
25
0 16.8 4.76
Matrix Factorization
Finding the [U] matrix
Using the Forward Elimination Procedure of Gauss Elimination
25 5 1 25 5 1
0 4.8 1.56 Row 2
Row 3 (16.8) 0 4.8 1.56
4.8
0 16.8 4.76
0
0 0.7
25 5 1
0 4.8 1.56
U
0
0 0. 7
Matrix Factorization
Finding the [V] matrix
Using the multipliers used during the Forward Elimination Procedure
a 21 64
From the first step 1 0 0 v21 2.56
of forward v 0
a11 25
elimination
21 1
v31 v32
1
a 31 144
v31 5.76
a 11 25
From the second 25 5 1
step of forward 0 4.8 1.56 a 32 16 .8
v32 3 .5
elimination 0 16.8 4.76
a 22 4.8
Matrix Factorization
1 0 0
V 2.56 1 0
5.76 3.5 1
Does V U A ?
1 0 0 25 5 1
V U 2.56 1 0 0 4.8 1.56
5.76 3.5 1 0
0 0.7
Matrix Factorization
Example: Solving simultaneous linear equations using VU factorization
Solve the following set of 25 5 1 a 1 106.8
linear equations using VU 64 8 1 a 177.2
2
Factorization 144 12 1 a 3 279.2
Using the procedure for finding the [V] and [U] matrices
1 0 0 25 5 1
A V U 2.56 1 0 0 4.8 1.56
5.76 3.5 1 0
0 0.7
Matrix Factorization
Example: Solving simultaneous linear equations using VU factorization
Complete the forward substitution to solve for Z
z1 106.8
z 2 177.2 2.56 z1
177.2 2.56(106.8) z1 106.8
96.2 z 96.21
Z 2
z 3 279.2 5.76 z1 3.5 z 2
279.2 5.76(106.8) 3.5(96.21) z3 0.735
0.735
Matrix Factorization
Example: Solving simultaneous linear equations using VU factorization
Set U X Z 25 5 1 a1 106.8
0 4.8 1.56 a - 96.21
2
0
0 0.7 a3 0.735
The 3 equations become
Solve for X 25a1 5a 2 a3 106.8
4.8a 2 1.56a3 96.21
0.7a3 0.735
Matrix Factorization
Example: Solving simultaneous linear equations using VU factorization
From the 3rd equation Substituting in a3 and using the
second equation
0.7a3 0.735
4.8a 2 1.56a3 96.21
0.735
a3 96.21 1.56a3
0.7 a2
4.8
1.050 - 96.21 1.561.050
- 4.8
19.70
Matrix Factorization
Example: Solving simultaneous linear equations using VU factorization
Substituting in a3 and a2 Hence the Solution Vector is:
using the first equation
25a1 5a2 a3 106.8 a1 0.2900
a1
106.8 5a 2 a 3 a 19.70
25 2
106.8 519.70 1.050 a3 1.050
25
0.2900
Gauss Method
An iterative method.
Basic Procedure:
-Algebraically solve each linear equation for xi
-Assume an initial guess solution array
-Solve for each xi and repeat
-Use absolute relative approximate error after each iteration
to check if error is within a prespecified tolerance.
Gauss Method
Why?
The Gauss Method allows the user to control round-off error.
Elimination methods such as Gaussian Elimination and VU
Factorization are prone to round-off error.
Also: If the physics of the problem is understood, a close initial
guess can be made, decreasing the number of iterations needed.
Gauss Method
Algorithm
A set of n equations and n unknowns:
If: the diagonal elements are
a11 x1 a12 x2 a13 x3 ... a1n xn b1 non-zero
a21 x1 a22 x2 a23 x3 ... a2n xn b2 Rewrite each equation solving
. .
. . for the corresponding unknown
. .
ex:
an1 x1 an 2 x2 an 3 x3 ... ann xn bn
First equation, solve for x1
Second equation, solve for x2
Gauss Method
Algorithm
Rewriting each equation
c a12 x 2 a13 x3 a1n x n From Equation 1
x1 1
a11
c2 a21 x1 a23 x3 a2 n xn
x2 From equation 2
a22
cn 1 an 1,1 x1 an 1, 2 x2 an 1,n 2 xn 2 an 1,n xn From equation n-1
xn 1
an 1,n 1
cn an1 x1 an 2 x2 an ,n 1 xn 1 From equation n
xn
ann
Gauss Method
Algorithm
General Form of each equation
n
a
n
c1 a1 j x j cn 1 n 1, j xj
j 1 j 1
j n 1
x1
j 1 xn 1
a11 an 1,n 1
n
c n a nj x j
n
c2 a2 j x j
j 1 j 1
j n
x2
j 2
xn
a 22 a nn
Gauss Method
Gauss Algorithm
General Form for any row ‘i’
n
ci aij x ( k 1) j
j 1
j i
x(k )i , i 1,2,, n.
aii
How or where can this equation be used?
Gauss Method
Solve for the unknowns
Assume an initial guess for [X] Use rewritten equations to solve for
each value of xi.
x1 Important: Remember to use the
most recent value of xi. Which
x means to apply values calculated to
2 the calculations remaining in the
current iteration.
xn -1
xn
Gauss Method
Calculate the Absolute Relative Approximate Error
x inew x iold
a i
new
100
xi
So when has the answer been found?
The iterations are stopped when the absolute relative
approximate error is less than a pre-specified tolerance for all
unknowns.
Gauss-Seidel Method
• Improved method
i 1 n
ci aij x j a x
( k 1)
(k )
ij j
j 1 j i 1
(k )
xi
aii
Gauss-Seidel Method: Pitfall
Diagonally dominant: The coefficient on the diagonal must be at least
equal to the sum of the other coefficients in that row and at least one row
with a diagonal coefficient greater than the sum of the other coefficients
in that row.
Which coefficient matrix is diagonally dominant?
2 5.81 34 124 34 56
A 45 43 1 [B] 23 53 5
123 16 1 96 34 129
Most physical systems do result in simultaneous linear equations that
have diagonally dominant coefficient matrices.
Gauss-Seidel Method: Example
Given the system of equations The coefficient matrix is:
12x 1 3x 2 - 5x 3 1
x 1 5x 2 3x 3 28 12 3 5
3x1 7x2 13x3 76
A 1 5 3
3 7 13
With an initial guess of
Will the solution converge using the
x1 1
Gauss-Seidel method?
x 0
2
x 3 1
Gauss-Seidel Method: Example
Checking if the coefficient matrix is diagonally dominant
a11 12 12 a12 a13 3 5 8
12 3 5
A 1 5 3 a22 5 5 a21 a23 1 3 4
3 7 13
a33 13 13 a31 a32 3 7 10
The inequalities are all true and at least one row is strictly greater than:
Therefore: The solution should converge using the Gauss-Seidel Method
Gauss-Siedel Method: Example
Rewriting each equation With an initial guess of
12 3 5 a1 1 x1 1
1 5 3 a 28 x 0
2 2
3 7 13 a3 76
x 3 1
1 3 x 2 5 x3 1 30 51
x1 x1 0.50000
12 12
28 x1 3x3 28 0.5 31
x2 x2 4.9000
5 5
76 3x1 7 x2
x3 76 30.50000 74.9000
13 x3 3.0923
13
Gauss-Siedel Method: Example
The absolute relative approximate error
0.50000 1.0000
a 1 100 67 .662 %
0.50000
4.9000 0
a 2
100 100 .00 %
4.9000
3.0923 1.0000
a 3
100 67 .662 %
3.0923
The maximum absolute relative error after the first iteration is 100%
Gauss-Siedel Method: Example
After Iteration #1
x1 0.5000
x 4.9000
2
x3 3.0923
Substituting the x values into the equations After Iteration #2
1 34.9000 53.0923 x1 0.14679
x1 0.14679 x 3.7153
12 2
x3 3.8118
28 0.14679 33.0923
x2 3.7153
5
76 30.14679 74.900
x3 3.8118
13
Gauss-Siedel Method: Example
Iteration #2 absolute relative approximate error
0.14679 0.50000
a 1 100 240 .62 %
0.14679
3.7153 4.9000
a 2 100 31 .887 %
3.7153
3.8118 3.0923
a 3 100 18 .876 %
3.8118
The maximum absolute relative error after the first iteration is 240.62%
This is much larger than the maximum absolute relative error obtained in
iteration #1. Is this a problem?
Gauss-Siedel Method: Example
Repeating more iterations, the following values are obtained
Iteration a1
a 1 a2
a 2
a3
a 3
1 0.50000 67.662 4.900 100.00 3.0923 67.662
2 0.14679 240.62 3.7153 31.887 3.8118 18.876
3 0.74275 80.23 3.1644 17.409 3.9708 4.0042
4 0.94675 21.547 3.0281 4.5012 3.9971 0.65798
5 0.99177 4.5394 3.0034 0.82240 4.0001 0.07499
6 0.99919 0.74260 3.0001 0.11000 4.0001 0.00000
x1 0.99919
The solution obtained
x2 3.0001
x3 4.0001
x1 1
is close to the exact solution of x 3
2
x 3 4
Comparison of different methods
• Gauss is slow convergent; however, more
stable.
• Gauss-Seidel might not be stable;
however, when stable, converges fast.
• Hotelling and bodewig uses an improved
iterative method for matrix inversion.
• In Relaxation techniques small corrections
will be applied to each element at each
iteration.
Transformations and Eigenvlues
• Many problems of the form: y Ax
• Can be written as: y Sx
• Where S is a diagonal Matrix.
• Where the primed vectors represent the original
vectors in the new space.
e d ij e j
i
j
e De
Eigenvalues and Eigenvectors
x D 1x
e De
y [DAD1 ]x Sx
DAD1 S
ADT D T S (Assuming Orthonorma transform
l ations)
Writing in component format :
n n
a
k 1
ik d jk d ji S jj d
k 1
jk ki S jj
n
(a
k 1
ik ki S jj ) d jk 0
Det aik ki S jj 0, j