Worksheet 4

In this worksheet we will provide you with some basics to understand the Gaussian Elimination.

We seek to solve the linear system given by \begin{equation} Ax = b, \end{equation} In general, to solve such systems we would like to transform the system to a more convenient form.
This is generally done by computing a factorization of the Matrix, LU,Cholesky, LDL, etc.

In the case of the Gaussian elimination, we transform $A$ to a upper triangular matrix. This process is equivalent to multiply the system by a lower triangular matrix \begin{equation} L^{-1}Ax = L^{-1}b, \end{equation} such that the resulting system is given by \begin{equation} Ux = L^{-1}b, \end{equation} where $U$ is an upper triangular matrix.

Gaussian elimination will perform the same operation using the augmented matrix (without the Matrix $L$).

Once the system is reduced to this form, we use a triangular solver (or backwards substitution to find $x$).
\begin{equation} x = U^{-1}( L^{-1}b). \end{equation}

Linear Solvers

Then we can divide the solution of $Ax = b$ in two stages:

  • Reduction to a triangular form (Gaussian Elimination);
  • Triangular solve (backward subtitution).

Gaussian Elimination

We present an example of Gaussian Elimination using a randomly generated $4\times 4$ matrix


In [20]:
n = 4
A = rand(n,n) + n*eye(n) # we add a identity to be sure that the matrix is invertible


Out[20]:
4x4 Array{Float64,2}:
 4.6126    0.00425747  0.976279  0.563531
 0.67515   4.292       0.567439  0.724985
 0.421752  0.344754    4.19598   0.116179
 0.331095  0.30775     0.613042  4.64796 

and we suppose a $5 \times 1$ right hand size $b$, given by


In [21]:
b = rand(n,1)


Out[21]:
4x1 Array{Float64,2}:
 0.750717
 0.668677
 0.448934
 0.573083

In order to perform the Gaussian elimination we build the augmented matrix


In [52]:
M = hcat(A,b) # or in equivalent fashion M = [A b]


Out[52]:
4x5 Array{Float64,2}:
 4.6126    0.00425747  0.976279  0.563531  0.750717
 0.67515   4.292       0.567439  0.724985  0.668677
 0.421752  0.344754    4.19598   0.116179  0.448934
 0.331095  0.30775     0.613042  4.64796   0.573083

We start the first step of the Gaussian elimination. We look for the pivot, i.e. the first non-zero in the first column


In [53]:
i = 1; #first column
find(M[i:end,i])


Out[53]:
4-element Array{Int64,1}:
 1
 2
 3
 4

The function find provides the indices of all the non zeros elements in the input vector. In this case all of them are non zero, so we don't need to perform a row swap. Then we can start introducing zeros in the first column, below the diagonal. In order to do this we perform a row operation, given by


In [54]:
M[2,:] = M[2,:] - M[2,1]/M[1,1]*M[1,:]
M


Out[54]:
4x5 Array{Float64,2}:
 4.6126    0.00425747  0.976279  0.563531  0.750717
 0.0       4.29138     0.424541  0.6425    0.558794
 0.421752  0.344754    4.19598   0.116179  0.448934
 0.331095  0.30775     0.613042  4.64796   0.573083

As you can see we introduced a zero in (2,1). We can repeat this operation with a for for the rest of the column


In [55]:
for j = 1+1:n
    M[j,:] = M[j,:] - M[j,1]/M[1,1]*M[1,:]
end
M


Out[55]:
4x5 Array{Float64,2}:
 4.6126  0.00425747  0.976279  0.563531   0.750717
 0.0     4.29138     0.424541  0.6425     0.558794
 0.0     0.344365    4.10671   0.0646523  0.380292
 0.0     0.307445    0.542965  4.60751    0.519197

You can see that we introduced zeros in the whole first column under the diagonal.
Now, we can start with the second row.
We check for the first non-zero in the second column, under the diagonal.


In [57]:
indnz = find(M[2:n,2])[1] # you pick the first index


Out[57]:
1

indz is the first non-zero index in the 2nd column starting from the 2nd row (i.e. is a local index) we need to transorm it back to a global index


In [58]:
indj = 1 + indnz


Out[58]:
2

If indj is different from i, this means that $M_{i,i} = 0$ and we need to change the pivot using a row swap


In [59]:
if indj != 2
    buffer = M[2,:]
    M[2,:] = M[indj,:]
    M[indj,:]  = buffer
end

Now we are sure that we have a non zero pivot, so we can continue to introduce zeros


In [60]:
for j = 3:n
    M[j,:] = M[j,:] - M[j,2]/M[2,2]*M[2,:]
end
M


Out[60]:
4x5 Array{Float64,2}:
 4.6126  0.00425747  0.976279  0.563531   0.750717
 0.0     4.29138     0.424541  0.6425     0.558794
 0.0     0.0         4.07265   0.0130944  0.335452
 0.0     0.0         0.512549  4.56148    0.479163

Now we can go to third column


In [61]:
indnz = find(M[3:n,3])[1] # you pick the first index
indj = 2 + indnz
if indj != 3
    buffer = M[3,:]
    M[3,:] = M[indj,:]
    M[indj,:] = buffer
end
for j = 4:n
    M[j,:] = M[j,:] - M[j,3]/M[3,3]*M[3,:]
end
M


Out[61]:
4x5 Array{Float64,2}:
 4.6126  0.00425747  0.976279  0.563531   0.750717
 0.0     4.29138     0.424541  0.6425     0.558794
 0.0     0.0         4.07265   0.0130944  0.335452
 0.0     0.0         0.0       4.55983    0.436946

Now we can extract $U$ and $L^{-1}b = b1$ from the augmented matrix


In [62]:
U = M[1:n,1:n]
b1 = M[1:n,n+1]


Out[62]:
4-element Array{Float64,1}:
 0.750717
 0.558794
 0.335452
 0.436946

Triangular Solve

From the Gaussian Eliminiation we have reduced the problem to solving \begin{equation} Ux = b, \end{equation} where U is an upper triangular matrix.

Now we need to solve for $x$, we create a vector $x$ of the same size as $b$


In [64]:
x = zeros(b)


Out[64]:
4x1 Array{Float64,2}:
 0.0
 0.0
 0.0
 0.0

we start by setting $x_{n,n} = b_{n}/U_{n,n}$


In [65]:
x[n] = b1[n]/U[n,n]


Out[65]:
0.09582498766476033

We solve the second to last equation that is given by \begin{equation} U_{n-1,n-1}x_{n-1} + U_{n-1,n}x_{n} = b_{n-1}, \end{equation} and we solve for $x_{n-1}$, which results in \begin{equation} x_{n-1} = \frac{b_{n-1}- U_{n-1,n}x_{n} }{U_{n-1,n-1}}, \end{equation} which can be translated to julia code as follows


In [67]:
x[n-1] = (b1[n-1] - U[n-1,n]*x[n])/U[n-1,n-1]


Out[67]:
0.08205887204535202

The third to last equation can be written as \begin{equation} U_{n-2,n-2}x_{n-2} + U_{n-2,n-1}x_{n-1} + U_{n-2,n}x_{n} = b_{n-2}, \end{equation} where we know $x_{n-1}$ and $x_{n}$ so we can solve $x_{n-2}$ using \begin{equation} x_{n-2} = \frac{ b_{n-2}- U_{n-2,n-1}x_{n-1} + U_{n-2,n}x_{n}}{U_{n-2,n-2}} \end{equation} where we can write $U_{n-2,n-1}x_{n-1} + U_{n-2,n}x_{n}$ as an inner product, which results in the following julia code


In [69]:
x[n-2] = (b1[n-2] - (U[n-2,n-2+1:n]*x[n-2+1:n])[1])/U[n-2,n-2]


Out[69]:
0.1077484335827115

The [1] in the dot product is needed to be consistent. The multiplication outputs a matrix with one element (but still a matrix), so you need to extract that one element out of the resulting matrix.
You can repeat the same for rest of the unknowns.


In [72]:
x[n-3] = (b1[n-3] - (U[n-3,n-3+1:n]*x[n-3+1:n])[1])/U[n-3,n-3]


Out[72]:
0.13357879992820912

Now we can check that our solution is the correct one, by checking the residual, which should be extremely small.


In [75]:
print("Residual = ",norm(A*x - b)/norm(b))


Residual = 1.3416725090382251e-16

Remark

Gaussian Elimination is analogous to compute the lu factorization of A. The lu factorization, finds two factors, one lower triangular $L$ and an upper triangular $U$, such that \begin{equation} A = LU. \end{equation}

If the LU factors are know, then we can compute the solution to $AX = b$ by a forward substition (i.e. we apply $L^{-1}$) and a backward substitution (i.e. we apply $U^{-1}$) which results in \begin{equation} x = U^{-1}( L^{-1}b). \end{equation}


In [ ]: