Open Educational Resources

Direct Methods: Naive Gauss Elimination

The Method

The naive Gauss elimination is a procedure in which the linear system of equations is manipulated such that the coefficients of the component x_i are eliminated from equation i+1 to equation n. It would then be possible to solve for x_n using the last equation, and then backward substitute to find the remaining components. The procedure is divided into two steps, forward elimination followed by backward substitution. We will illustrate the process in the following example. Consider the system of equations Ax=b defined as:

    \[\begin{split} 2x_1+2x_2+2x_3&=6\\ 2x_1+3x_2+4x_3&=4\\ -x_1+2x_2+3x_3&=8\\ \end{split} \]

Forward elimination: First, we will use the first equation to eliminate the coefficients of x_1 in equation 2 and 3. By multiplying the first equation by -1 and adding it to the second equation we get:

    \[\begin{split} 2x_1+2x_2+2x_3&=6\\ (2-2)x_1+(3-2)x_2+(4-2)x_3&=(4-6)\\ -x_1+2x_2+3x_3&=8\\ \end{split} \]

which becomes:

    \[\begin{split} 2x_1+2x_2+2x_3&=6\\ 0+1x_2+2x_3&=-2\\ -x_1+2x_2+3x_3&=8\\ \end{split} \]

By multiplying the first equation by \frac{1}{2} and adding it to the third equation we get:

    \[\begin{split} 2x_1+2x_2+2x_3&=6\\ 0+1x_2+2x_3&=-2\\ 0+3x_2+4x_3&=11\\ \end{split} \]

We are now ready to eliminate the coefficient of x_2 from equation 3 by multiplying equation 2 by -3 and adding it to equation 3:

    \[\begin{split} 2x_1+2x_2+2x_3&=6\\ 0+1x_2+2x_3&=-2\\ 0+0-2x_3&=17\\ \end{split} \]

Backward substitution: Solving the above system becomes very simple. x_3 is straightforward, it is equal to x_3=\frac{17}{-2}=-8.5. x_2 can be obtained from the second equation as: x_2=-2-2x_3=-2-2(-8.5)=15 and finally x_1=\frac{6-2(-8.5)-2(15)}{2}=-3.5. Therefore, the solution is:

    \[ x=\left(\begin{array}{c}-3.5\\15\\-8.5\end{array}\right) \]

Of course, as expected, if we multiply the original A by x we should get the matrix b.

The Naive Gauss elimination method can be implemented in matrix form of dimensions n\times (n+1) as follows:

    \[ \left(\begin{array}{ccc|c}2&2&2&6\\2&3&4&4\\-1&2&3&8\end{array}\right)\Rightarrow \left(\begin{array}{ccc|c}2&2&2&6\\0&1&2&-2\\-1&2&3&8\end{array}\right) \Rightarrow \left(\begin{array}{ccc|c}2&2&2&6\\0&1&2&-2\\0&3&4&11\end{array}\right) \Rightarrow \left(\begin{array}{ccc|c}2&2&2&6\\0&1&2&-2\\0&0&-2&17\end{array}\right) \]

The following code implements the above calculations in Mathematica

View Mathematica Code
A = {{2, 2, 2}, {2, 3, 4}, {-1, 2, 3}};
b = {6, 4, 8};
Aa = {{2, 2, 2, b[[1]]}, {2, 3, 4, b[[2]]}, {-1, 2, 3, b[[3]]}};
Aa[[2]] = Aa[[2]] + -1*Aa[[1]]
Aa[[3]] = Aa[[3]] + 1/2*Aa[[1]]
Aa[[3]] = Aa[[3]] - 3 Aa[[2]]
Aa // MatrixForm
x3 = Aa[[3, 4]]/Aa[[3, 3]]
x2 = (Aa[[2, 4]] - Aa[[2, 3]]*x3)/Aa[[2, 2]]
x1 = (Aa[[1, 4]] - Aa[[1, 2]]*x2 - Aa[[1, 3]]*x3)/Aa[[1, 1]]
View Python Code
import numpy as np
A = [[2, 2, 2], [2, 3, 4], [-1, 2, 3]]
b = [6, 4, 8]
Aa = np.array([[2, 2, 2, b[0]], [2, 3, 4, b[1]], [-1, 2, 3, b[2]]])
Aa[1] = Aa[1] + -1*Aa[0]
print("Aa[1]:",Aa[1])
Aa[2] = Aa[2] + 1/2*Aa[0]
print("Aa[2]:",Aa[2])
Aa[2] = Aa[2] - 3*Aa[1]
print("Aa[2]:",Aa[2])
print("Aa:",Aa)
x3 = Aa[2, 3]/Aa[2, 2]
print("x3:",x3)
x2 = (Aa[1, 3] - Aa[1, 2]*x3)/Aa[1, 1]
print("x2:",x2)
x1 = (Aa[0, 3] - Aa[0, 1]*x2 - Aa[0, 2]*x3)/Aa[0, 0]
print("x1:",x1)

Programming the Method

The following are the steps to program the Naive Gauss Elimination method. Assuming an n number of equations in n unknowns of the form Ax=b:

  1. Form the combined matrix G=A|b
  2. Forward elimination: Use the pivot elements G_{kk} on the row k as “Pivot” elements. Use the “Pivot” elements to eliminate the components G_{ik} with i from k+1 to n. This is done by iterating i from k+1 to n and for each row i, doing the operation row_i=row_i-row_k\frac{G_{ik}}{G_{kk}}.
  3. Backward substitution: The element x_i with i running backwards from n to 1 can be found using the equation:

        \[ x_i=\frac{G_{i(n+1)}-\sum_{l=i+1}^nG_{il}x_l}{G_{ii}} \]

The following is the code that relies on the “Sum” built-in function in Mathematica.

View Mathematica Code
(*Procedure*)
GaussElimination[A_, b_] :=
 (n = Length[A];
  G = Transpose[Insert[Transpose[A], b, n + 1]];
  For[k = 1, k <= n, k++, 
   Do[G[[i]] = G[[i]] - G[[i, k]]/G[[k, k]]*G[[k]], {i, k + 1, n}]];
  x = Table[0, {i, 1, n}];
  For[i = n, i >= 1, i = i - 1, 
   x[[i]] = (G[[i, n + 1]] - Sum[G[[i, l]]*x[[l]], {l, i + 1, n}])/G[[i, i]]];
  x)
(*An example*)
A = {{1, 2, 3, 5}, {2, 0, 1, 4}, {1, 2, 2, 5}, {4, 3, 2, 2}};
b = {-4, 8, 0, 10};
(*Applying the procedure to the example*)
N[GaussElimination[A, b]]
GaussElimination[A, b]
View Python Code
import numpy as np
# Procedure
def GaussElimination(A,b):
  n = len(A)
  G = (np.vstack([A.astype(np.float).T, b.astype(np.float)])).T
  for k in range(n): 
    for i in range(k + 1, n):
      G[i] = G[i] - G[i][k]/G[k][k]*G[k]
  x = np.zeros(n)
  for i in range(n-1,-1,-1): 
    x[i] = (G[i][n] - sum([G[i][l]*x[l] for l in range(i + 1, n)]))/G[i][i]
  return x

# An example
A = np.array([[1, 2, 3, 5], [2, 0, 1, 4], [1, 2, 2, 5], [4, 3, 2, 2]])
b = np.array([-4, 8, 0, 10])
# Applying the procedure to the example
GaussElimination(A, b)

The following is the code that does not rely on the built-in function in Mathematica.

View Mathematica Code
(*Procedure*)
GaussElimination[A_, b_] := (n = Length[A];
  G = Transpose[Insert[Transpose[A], b, n + 1]];
  For[k = 1, k <= n, k++, 
   Do[G[[i]] = G[[i]] - G[[i, k]]/G[[k, k]]*G[[k]], {i, k + 1, n}]];
  x = Table[0, {i, 1, n}];
  
  For[i = n, i >= 1, 
   i = i - 1, (BB = 0; Do[BB = G[[i, l]]*x[[l]] + BB, {l, i + 1, n}]; 
    x[[i]] = (G[[i, n + 1]] - BB)/G[[i, i]])];
  x)
(*An example*)
A = {{1, 2, 3, 5}, {2, 0, 1, 4}, {1, 2, 2, 5}, {4, 3, 2, 2}};
b = {-4, 8, 0, 10};
(*Applying the procedure to the example*)
N[GaussElimination[A, b]]
GaussElimination[A, b]
View Python Code
import numpy as np
# Procedure
def GaussElimination(A,b):
  n = len(A)
  G = (np.vstack([A.astype(np.float).T, b.astype(np.float)])).T
  for k in range(n): 
    for i in range(k + 1, n):
      G[i] = G[i] - G[i][k]/G[k][k]*G[k]
  x = np.zeros(n)
  for i in range(n-1,-1,-1):
    BB = 0
    for l in range(i+1, n):
      BB = G[i][l]*x[l] + BB
    x[i] = (G[i][n] - BB)/G[i][i]
  return x

# An example
A = np.array([[1, 2, 3, 5], [2, 0, 1, 4], [1, 2, 2, 5], [4, 3, 2, 2]])
b = np.array([-4, 8, 0, 10])
# Applying the procedure to the example
GaussElimination(A, b)

The following link provides the MATLAB code for implementing the Naive Gauss Elimination Method.

Why Naive

The method described above is considered “naive” for two reasons. The first reason is that it fails to find a solution if a pivot element (one of the diagonal elements) is found to be zero. For example, consider the linear system:

    \[ \left(\begin{matrix}1&0&0\\0&0&1\\0&1&0\end{matrix}\right)\left(\begin{matrix}x_1\\x_2\\x_3\end{matrix}\right)=\left(\begin{matrix}1\\1\\1\end{matrix}\right) \]

This is really simple and the solution is just x_1=x_2=x_3=1. However, the code fails to find a solution because when the second row is used for elimination, the pivot element A_{22} is equal to zero, and so, an error is encountered. Try the above systems out and see the error that Mathematica spits out. This problem can be solved by simply exchanging rows 2 and 3. The same issue will be encountered for the following system:

    \[ \left(\begin{matrix}1&1&0\\0&0&1\\1&0&0\end{matrix}\right)\left(\begin{matrix}x_1\\x_2\\x_3\end{matrix}\right)=\left(\begin{matrix}1\\1\\1\end{matrix}\right) \]

In this case, one solution would be to rearrange the equations such that equation 3 would be equation 1. In general, the naive Gauss elimination code presented above can be augmented so that the pivot elements are chosen to be the largest elements in the column below. This augmentation is termed: “Partial Pivoting”. Partial Pivoting is performed by exchanging the rows. In this case, the order of the variables x_i stays the same. In rare cases, there could be a need to exchange the columns as well; a procedure termed “Complete Pivoting”, but this adds the complication of re-arranging the variables.

The second reason why it is called naive is that the method essentially produces an upper triangular matrix that can be easily solved. It is important to note that if the system is already a lower triangular matrix, then there is no need to employ the forward elimination procedure, because we can simply use forward substitution to find x_1, x_2, \cdots,x_n in this order.
Consider the following system:

    \[ \left(\begin{matrix}1&0&0\\1&1&0\\1&1&1\end{matrix}\right)\left(\begin{matrix}x_1\\x_2\\x_3\end{matrix}\right)=\left(\begin{matrix}1\\0\\1\end{matrix}\right) \]

It is straightforward that x_1=1, x_2=0-x_1=-1, and x_3=1-x_1-x_2=1-1+1=1. However, employing the naive Gauss elimination would conduct the following steps:

    \[ \left(\begin{array}{ccc|c}1&0&0&1\\1&1&0&0\\1&1&1&1\end{array}\right)\Rightarrow \left(\begin{array}{ccc|c}1&0&0&1\\0&1&0&-1\\1&1&1&1\end{array}\right) \Rightarrow \left(\begin{array}{ccc|c}1&0&0&1\\0&1&0&-1\\0&1&1&0\end{array}\right) \Rightarrow \left(\begin{array}{ccc|c}1&0&0&1\\0&1&0&-1\\0&0&1&1\end{array}\right) \]

Leave a Reply

Your email address will not be published. Required fields are marked *