Gaussian elimination

We create a 4-by-4 linear system with the matrix

A = [
     2    0    4    3 
    -4    5   -7  -10 
     1   15    2   -4.5
    -2    0    2  -13
    ];

and with the right-hand side

b = [ 4, 9, 29, 40 ];

We define an augmented matrix by tacking \(b\) on the end as a new column.

S = [A b]
4×5 Array{Float64,2}:
  2.0   0.0   4.0    3.0   4.0
 -4.0   5.0  -7.0  -10.0   9.0
  1.0  15.0   2.0   -4.5  29.0
 -2.0   0.0   2.0  -13.0  40.0

The goal is to introduce zeros into the lower triangle of this matrix. By using only elementary row operations, we ensure that the matrix \(S\) always represents a linear system that is equivalent to the original. We proceed from left to right and top to bottom. The first step is to put a zero in the (2,1) location using a multiple of row 1:

@show mult21 = S[2,1]/S[1,1];
S[2,:] -= mult21*S[1,:];   # -= means "subtract the following from"
S
mult21 = S[2, 1] / S[1, 1] = -2.0
4×5 Array{Float64,2}:
  2.0   0.0  4.0    3.0   4.0
  0.0   5.0  1.0   -4.0  17.0
  1.0  15.0  2.0   -4.5  29.0
 -2.0   0.0  2.0  -13.0  40.0

We repeat the process for the (3,1) and (4,1) entries.

@show mult31 = S[3,1]/S[1,1];
S[3,:] -= mult31*S[1,:];
@show mult41 = S[4,1]/S[1,1];
S[4,:] -= mult41*S[1,:];
S
mult31 = S[3, 1] / S[1, 1] = 0.5
mult41 = S[4, 1] / S[1, 1] = -1.0
4×5 Array{Float64,2}:
 2.0   0.0  4.0    3.0   4.0
 0.0   5.0  1.0   -4.0  17.0
 0.0  15.0  0.0   -6.0  27.0
 0.0   0.0  6.0  -10.0  44.0

The first column has the zero structure we want. To avoid interfering with that, we no longer add multiples of row 1 to anything. Instead, to handle column 2, we use multiples of row 2. We’ll also exploit the highly repetitive nature of the operations to write them as a loop.

for i = 3:4
    mult = S[i,2]/S[2,2]
    S[i,:] -= mult*S[2,:]
end
S
4×5 Array{Float64,2}:
 2.0  0.0   4.0    3.0    4.0
 0.0  5.0   1.0   -4.0   17.0
 0.0  0.0  -3.0    6.0  -24.0
 0.0  0.0   6.0  -10.0   44.0

We finish out the triangularization with a zero in the (4,3) place. It’s a little silly to use a loop for just one iteration, but the point is to establish a pattern.

for i = 4
    mult = S[i,3]/S[3,3]
    S[i,:] -= mult*S[3,:]
end
S
4×5 Array{Float64,2}:
 2.0  0.0   4.0   3.0    4.0
 0.0  5.0   1.0  -4.0   17.0
 0.0  0.0  -3.0   6.0  -24.0
 0.0  0.0   0.0   2.0   -4.0

Recall that \(S\) is an augmented matrix: it represents the system \(Ux=z\), where

U = S[:,1:4]
4×4 Array{Float64,2}:
 2.0  0.0   4.0   3.0
 0.0  5.0   1.0  -4.0
 0.0  0.0  -3.0   6.0
 0.0  0.0   0.0   2.0
z = S[:,5]
4-element Array{Float64,1}:
   4.0
  17.0
 -24.0
  -4.0

The solutions to this system are identical to those of the original system, but this one can be solved by backward substitution.

using FundamentalsNumericalComputation
x = FNC.backsub(U,z)
4-element Array{Float64,1}:
 -3.0
  1.0
  4.0
 -2.0
b - A*x
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0