Systems of Linear Equations

1
Systems of Linear Equations
2
Systems of Linear Equations

• Vectors
• Matrices
• Special matrices
• Vector vector scalar product
• Matrix multiplied by vector
• Linear system of equations
• Matrix multiplication

3
Systems of Linear Equations

• Is there a solution?
• If there is a solution, is there only one?
• How can we express the solution with the inverse
  of the coefficient matrix?
• Is it usually a good idea to get the solution via
  the inverse?

4

• The determinant is a number that can be
  calculated knowing the elements of the
  (cofficient) matrix.
• The matrix inverse doesnt exist if the
  determinant is zero.
• A condition for the inverse to exist is that the
  determinant is not zero.
• Determinant zero coefficient matrix is singular
  solution does not exist system is contradictory
  OR undetermined matrix rank is less than number
  of equations (or just division by zero, overflow
  etc.) all mean the same

5
Solution of Systems of Linear Equations

• Ax b
• solution
• X A-1b (A-1 is the inverse matrix)
• or, Ix A-1b (I is identity matrix)

6
Types of matrix problems
7
Example 1 - Unique solution.

• Consider the system
• 2x 3y 8
• 5x - 6y 30
• Geometrically this corresponds to two lines
  crossing in one point.

8
Example 2 - No unique solution.

• Consider the system
• 2x 3y 8
• 4x 6y 16
• Geometrically this corresponds to two coincident
  lines. Infinite number of solutions.

9
Example 3 - No solution.

• Consider the system
• 2x 3y 8
• 4x 6y 30
• Geometrically this corresponds to two parallel
  lines never crossing.

10
Example 4 Ill-conditioned matrix

• Consider the system
• 2x 3y 8
• 4x 7y 15
• Geometrically this corresponds to two almost
  parallel lines barely crossing. Sensitive to
  round-off error.

11
Existence and Uniqueness Statement with Rank

• The rank of a matrix tells us how many
  independent rows (or how many independent
  columns) the coefficient matrix has -- it may be
  less than the total number of equations.
• If rank is less than the number of equations, two
  cases can happen
• Consistent (if rank of coeff matrix rank of
  augmented matrix) Infinite number of solutions
• Non-Consistent (if rank of coeff matrix lt rank of
  augmented matrix) No solution

12
Definition Augmented matrix
13
Example 1

• Consider the system
• 2x 3y 8
• 4x 6y 16
• Thus the determinant is 26 - 34 0.
• And indeed the second equation is 2 times the
  first.
• The solution can be written as x 4 - 3/2 y for
  any y.
• Geometrically this corresponds to two coincident
  lines.
• In terms of rank the coefficient matrix has rank
  1, the augmented matrix has rank 1

14
Example 2

• Consider the system
• 2x 3y 8
• 4x 6y 17
• Thus the determinant is 26 - 34 0.
• There is no solution.
• Geometrically this corresponds to two parallel
  lines never crossing.
• In terms of rank the coefficient matrix has rank
  1, the augmented matrix has rank 2

15
Direct Methods

• Gauss-Jordan elimination
• Gaussian elimination
• LU decomposition / factorization
• Thomas algorithm

Iterative Methods

• Gauss-Seidel

16
Gaussian-Jordan

• Step 1. Forward elimination
• Choose a pivot element on the main diagonal
  (start at top left corner)
• Eliminate (set to zero) all elements above and
  below the main diagonal
• Move to the next column
• Repeat until matrix is an Identity matrix

17
Gauss-Jordan Elimination

• Reduce an augmented matrix to the identity to
  solve a system of linear equations.
• Example

18
Gauss-Jordan Elimination
19
Gauss-Jordan Elimination
20
Gauss-Jordan Elimination
21
Gauss-Jordan Elimination
22
Gauss-Jordan Elimination
23
Gauss-Jordan Elimination
24
Summary

• Gauss-Jordan process
• loop over rows (i)
• scale row such that diagonal is 1.0
• subtract a multiple of row i from every other row
  in order to fill the rest of the column with
  zeroes.
• move to next row
• Final solution appears at the end of the the
  augmented matrix.

25
Sub GaussJordan(A() As Double)
'Augmented matrix Dim PivElt As Double, TarElt As
Double Dim n As Integer
'number of equations Dim PivRow As Integer,
TarRow As Integer 'pivot row target row Dim i
As Integer, j As Integer n UBound(A, 1) For
PivRow 1 To n 'process
every row PivElt A(PivRow, PivRow)
'choose pivot element If PivElt 0 Then
MsgBox ("Zero pivot element encountered")
Stop End If For j 1 To n 1
A(PivRow, j) A(PivRow, j) / PivElt 'divide
whole row Next For TarRow 1 To n
'now replace all other rows If
Not (TarRow PivRow) Then TarElt
A(TarRow, PivRow) For j 1 To n 1
A(TarRow, j) A(TarRow, j) - A(PivRow, j)
TarElt Next j End If Next TarRow
Next PivRow End Sub
26
Gaussian Elimination

• Pivoting
• If the main diagonal pivot element is small or
  zero, interchange rows (from below) to maximize
  the pivot element

• The determinant
• After forward elimination, the determinant is the
  product of all the main diagonal elements of the
  upper triangular matrix. det(A) u11?u22?u33
  ??? ?unn

27
Gaus. Elim. forward elimination

• Eliminate elements in column 1
• Example

28
Gaus. Elim. forward elimination

• Eliminate a21 by calculating a ratio R
  a21/a11
• Set row 2 row 2 row 1

29
Gaus. Elim. forward elimination

• Continue until all elements below main diagonal
  0
• Matrix is now upper triangular

30
Gaus. Elim. backward substitution
31
Gaussian elimination

• Step 1. Forward elimination
• Choose a pivot element on the main diagonal
  (start at top left corner)
• Eliminate (set to zero) all elements below the
  main diagonal
• Move to the next column
• Repeat until matrix is upper triangular

• Step 2. Backward substitution
• Solve the last row explicitly for the last
  unknown
• Solve explicitly for the next row up
• Repeat until all the rows are solved

32
LU Factorization

• LU Factorization is a form of Gaussian
  Elimination
• The A matrix (in the problem Axb) is factored
  into the product of two matrices L and U (ALU)
• L is a lower triangular matrix and U is an upper
  triangular matrix

33
LU Factorization

• How do we find the entries in L and U?
• there are nine unknowns (3 Ls and 6 Us)
• there are nine equations (each setting a row of L
  times a column of U equal to one of the A values)

34
LU Factorization
35
Pivoting

• What would happen if a11 were zero?
• We would get a divide by zero error when building
  L.

• Algorithms for LU factorization check for this
  and swap the order of the rows in A (pivoting) to
  avoid this problem.

36
LU Factorization

• Once L and U have been determined we have
• LUx b

• How do we solve for x?
• Set UxL-1b y
• Therefore Ly b
• Solve Lyb for y
• Solve Ux y for x

37
Forward Substitution

• We solve Lyb using a process called forward
  substitution.

38
Back Substitution

• We solve Uxy using a process called back
  substitution.

39
Back Substitution
40
LU Decomposition and Backsubstitution
Basic Call LUDecompose(A, Indx) Call
LUSubstitute(A, Indx, b, x)
41
Sub LUDecompose(A() As Double, Indx() As
Integer) ' input a(n,n) ' output a(n,n) and
indx(n) Const Tiny As Double 1E-16 Dim n As
Integer n UBound(A, 1) Dim i As Integer, j As
Integer, k As Integer, m As Integer Dim dum As
Double For k 1 To n - 1 m k For i k
1 To n If (Abs(A(i, k)) gt Abs(A(m, k)))
Then m i Next i Indx(k) m dum A(m,
k) If m gt k Then A(m, k) A(k, k)
A(k, k) dum End If If Abs(dum) gt Tiny
Then dum 1 / dum Else MsgBox
("Singular coefficient matrix") Stop End
If Next k End Sub
42
Sub LUSubstitute(A() As Double, Indx() As
Integer, _ b() As Double, x() As
Double) ' input A(n,n) containing the LU
decomposition ' input Indx(n) index vector '
input b(n) right hand side ' output x(n)
solution Dim n As Integer n UBound(A, 1) Dim i
As Integer, j As Integer, k As Integer Dim dum As
Double For i 1 To n End sub
43
Direct vs Iterative

• Gauss elimination, Gauss-Jordan, LU-decomposition
  are DIRECT methods
• We can predict the number of multiplications and
  additions if the number of equations (n) is known
  
• There are other direct methods for special
  matrices. One of them is important for us It is
  the Thomas algorithm for Tridiagonal equations.

44
Tri-diagonal matrix

• Occurs often in reservoir simulation
• Thomas algorithm is efficient
• Uses LU factorization

45
Tri-diagonal matrix
X1 . . . . . . . . . . . . . xn
d1 . . . . . . . . . . . . . dn
46
Thomas algorithm save a,b,c vectors, not aii
elements
X1 . . . . . . . . . . . . . xn
d1 . . . . . . . . . . . . . dn
ai
bi
ci
47
Thomas algorithm
Sub Thomas (a() As Double, b() As Double, c() As
Double, d() As Double, x() As Double) 'Tridiagonal
system of equations Dim n As Integer, i As
Integer n UBound(b) ReDim w(n) As Double, g(n)
As Double w(1) b(1) g(1) d(1) / w(1)
For i 2 To n w(i) b(i) - a(i) c(i
- 1) / w(i - 1) g(i) (d(i) - a(i)
g(i - 1)) / w(i) Next i x(n) g(n) For i
n - 1 To 1 Step -1 x(i) g(i) - c(i) x(i
1) / w(i) Next i End Sub
48
Direct vs Iterative

• Iterative methods are generalization of the
  direct iteration we learned for one variable
  xg(x) equations.
• Iterative methods are used for very large but
  sparse systems (reservoir simulation).
• You can not predict the number of operations for
  an indirect method and you need a starting guess
  and a convergence criterion.
• The method we are going to learn is the
  Gauss-Seidel method.

49
Iterative methods

• Gauss Seidel

50
Gauss-Seidel

• Suppose we are trying to solve a system of three
  equations

51
Gauss-Seidel

• These equations can be rewritten as

52
Gauss-Seidel Iteration

• Gauss-Seidel starts from an initial guess of the
  xs and updates that guess in an attempt to
  converge to a solution

53
Diagonally Dominant Matrices

• Diagonally dominant matrices
• The absolute value of the diagonal element is
  greater than the sum of the absolute values of
  all the other elements in a given row.
• Identity matrix is diagonally dominant
• Diagonally Dominant?
• YES NO

54
Diagonal Dominance

• Gauss-Seidel will only work if/only if the matrix
  is diagonally dominant.