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}
Then we can divide the solution of $Ax = b$ in two stages:
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]:
and we suppose a $5 \times 1$ right hand size $b$, given by
In [21]:
b = rand(n,1)
Out[21]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
In [64]:
x = zeros(b)
Out[64]:
we start by setting $x_{n,n} = b_{n}/U_{n,n}$
In [65]:
x[n] = b1[n]/U[n,n]
Out[65]:
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]:
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]:
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]:
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))
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 [ ]: