Gaussian Elimination to Solve Linear Equations

The Gaussian Elimination Method is a widely used technique for solving systems of linear equations, where multiple equations with unknown variables are solved simultaneously to find the values of the unknowns. This method has practical applications in real life, such as in traffic flow analysis, where it helps solve systems representing traffic movement at intersections, optimizing the flow of traffic and reducing congestion.

In LU decomposition, Gaussian elimination is used to convert the matrix A into an upper triangular matrix (U).

Categories of Linear Equation Systems:

• Consistent Independent System: Has exactly one solution.
• Consistent Dependent System: Has infinite solutions.
• Inconsistent System: Has no solution.

Gaussian Elimination Method

Gaussian elimination is a row reduction algorithm for solving linear systems. It involves a series of operations on the augmented matrix (which includes both coefficients and constants) to simplify it into a row echelon form or reduced row echelon form. This method can also help in determining the rank, determinant and inverse of matrices. Gaussian elimination is a method for solving systems of equations in matrix form.

Let's for example we have system of linear equations:

a_1x+b_1y+c_1z=d_1
a_2x+b_2y+c_2z=d_2
a_3x+b_3y+c_3z=d_3
Matrix form of this system of linear equations: \begin{bmatrix}a_1 & b_1 & c_1 & |&d_1 \\a_2 & b_2 & c_2 & |&d_2 \\a_3 & b_3 & c_3 &| &d_3\end{bmatrix}

Goal:

Turn matrix into row-echelon form like, \begin{bmatrix}0 & a & b & |&d \\0 & 1 & c & |&e \\0 & 0 & 1 &| &f\end{bmatrix} .

Once in this form, we can say that 𝑧 = 𝑓 and use back substitution to solve for y and x.

Elementary Row Operations

1. Interchanging Rows: Swap two rows.
2. Multiplying a Row by a Scalar: Multiply all elements of a row by a non-zero number.
3. Adding a Scalar Multiple of One Row to Another: Add or subtract a multiple of one row to/from another.

These operations simplify solving systems of linear equations without changing the solution set.

Follow these steps:

1. Get a 1 in the first column, first row
2. Use the 1 to get 0’s in the remainder of the first column
3. Get a 1 in the second column, second row
4. Use the 1 to get 0’s in the remainder of the second column
5. Get a 1 in the third column, third row.

The resulting matrix is in row echelon form. A matrix is in reduced row-echelon form if all leading coefficients are 1 and the columns containing these leading coefficients have zeros elsewhere. This final form is unique, regardless of the sequence of row operations used. An example below illustrates this concept.

Example

To solve the system of equations:

1. x+y=3
2. 3x-2y=4

Step 1: Write the Coefficient Matrix

The coefficients of the unknowns (x and y) are written as a matrix:

\begin{bmatrix} 1 & 1 \\ 3 & -2 \end{bmatrix}

Step 2: Write the Augmented Matrix

Include the constants from the right-hand sides of the equations as an additional column:

\begin{bmatrix} 1 & 1 & |&3 \\ 3 & -2 & |&4 \end{bmatrix}

Step 3: Perform Row Operations

The goal is to simplify the augmented matrix to solve for x and y.

Subtract 3 \times r_1 from r_2 (to eliminate the x-term in the second row): r_2→r_2−3r_1

Calculation:\begin{bmatrix} 1 & 1 & 3 \\ 3 & -2 & 4 \end{bmatrix} \to \begin{bmatrix} 1 & 1 & 3 \\ 0 & -5 & -5 \end{bmatrix}

Simplify the second row:
Divide r_2​ by −5 to solve for y: \to \frac{r_2}{-5}

Result: \begin{bmatrix} 1 & 1 & 3 \\ 0 & 1 & 1 \end{bmatrix}

Step 4: Back-Substitute

Use the second row (y=1) to substitute into the first row (x + y = 3): x + 1 = 3  ⟹  x= 2

Solution

The solution to the system is: (x, y )=(2, 1)

Other algorithm in Gaussian elimination to solve Linear Equations

1. Partial Pivoting: Select the pivot element with the largest absolute value from the column. Swap rows to move this element into the pivot position. This improves numerical stability by reducing rounding errors.
2. Elimination: For each row below the pivot, subtract a multiple of the pivot row to eliminate the entries below the pivot.

Implementation of Gaussian Elimination to solve Linear Equations

1. Back Substitution

We use back-substitution to find the values of the unknowns after obtaining the row echelon form. The steps are:

1. Start with the last row of the row echelon form.
2. Solve for the last variable by dividing the constant term by its coefficient.
3. Substitute the known value of the last variable into the previous rows.
4. Solve for the next variable by isolating it and using the known values of the other variables.
5. Repeat the process for all rows, moving upwards to the first row.
6. Once all variables are solved, the system is complete and the solution is found.

Below is the implementation of the above algorithm.  

#include <bits/stdc++.h> using namespace std; #define N 3 int forwardElim(double mat[N][N + 1]); void backSub(double mat[N][N + 1]); void gaussianElimination(double mat[N][N + 1]) {  int singular_flag = forwardElim(mat);  if (singular_flag != -1)  {  printf("Singular Matrix.\n");  if (mat[singular_flag][N])  printf("Inconsistent System.");  else  printf("May have infinitely many solutions.");  return;  }  backSub(mat); } void swap_row(double mat[N][N + 1], int i, int j) {  for (int k = 0; k <= N; k++)  {  double temp = mat[i][k];  mat[i][k] = mat[j][k];  mat[j][k] = temp;  } } void print(double mat[N][N + 1]) {  for (int i = 0; i < N; i++, printf("\n"))  for (int j = 0; j <= N; j++)  printf("%lf ", mat[i][j]);  printf("\n"); } int forwardElim(double mat[N][N + 1]) {  for (int k = 0; k < N; k++)  {  int i_max = k;  int v_max = mat[i_max][k];  for (int i = k + 1; i < N; i++)  if (abs(mat[i][k]) > v_max)  v_max = mat[i][k], i_max = i;  if (!mat[k][i_max])  return k;  if (i_max != k)  swap_row(mat, k, i_max);  for (int i = k + 1; i < N; i++)  {  double f = mat[i][k] / mat[k][k];  for (int j = k + 1; j <= N; j++)  mat[i][j] -= mat[k][j] * f;  mat[i][k] = 0;  }  }  return -1; } void backSub(double mat[N][N + 1]) {  double x[N];  for (int i = N - 1; i >= 0; i--)  {  x[i] = mat[i][N];  for (int j = i + 1; j < N; j++)  {  x[i] -= mat[i][j] * x[j];  }  x[i] = x[i] / mat[i][i];  }  printf("\nSolution for the system:\n");  for (int i = 0; i < N; i++)  printf("%lf\n", x[i]); } int main() {  double mat[N][N + 1] = {{3.0, 2.0, -4.0, 3.0}, {2.0, 3.0, 3.0, 15.0}, {5.0, -3, 1.0, 14.0}};  gaussianElimination(mat);  return 0; } 

#include <bits/stdc++.h> using namespace std; #define N 3 int forwardElim(double mat[N][N + 1]); void backSub(double mat[N][N + 1]); void gaussianElimination(double mat[N][N + 1]) {  int singular_flag = forwardElim(mat);  if (singular_flag != -1) {  printf("Singular Matrix.\n");  if (mat[singular_flag][N])  printf("Inconsistent System.");  else  printf("May have infinitely many solutions.");  return;  }  backSub(mat); } void swap_row(double mat[N][N + 1], int i, int j) {  for (int k = 0; k <= N; k++) {  double temp = mat[i][k];  mat[i][k] = mat[j][k];  mat[j][k] = temp;  } } void print(double mat[N][N + 1]) {  for (int i = 0; i < N; i++, printf("\n"))  for (int j = 0; j <= N; j++)  printf("%lf ", mat[i][j]);  printf("\n"); } int forwardElim(double mat[N][N + 1]) {  for (int k = 0; k < N; k++) {  int i_max = k;  int v_max = mat[i_max][k];  for (int i = k + 1; i < N; i++)  if (abs(mat[i][k]) > v_max)  v_max = mat[i][k], i_max = i;  if (!mat[k][i_max])  return k;  if (i_max != k)  swap_row(mat, k, i_max);  for (int i = k + 1; i < N; i++) {  double f = mat[i][k] / mat[k][k];  for (int j = k + 1; j <= N; j++)  mat[i][j] -= mat[k][j] * f;  mat[i][k] = 0;  }  }  return -1; } void backSub(double mat[N][N + 1]) {  double x[N];  for (int i = N - 1; i >= 0; i--) {  x[i] = mat[i][N];  for (int j = i + 1; j < N; j++) {  x[i] -= mat[i][j] * x[j];  }  x[i] = x[i] / mat[i][i];  }  printf("\nSolution for the system:\n");  for (int i = 0; i < N; i++)  printf("%lf\n", x[i]); } int main() {  double mat[N][N + 1] = { { 3.0, 2.0, -4.0, 3.0 },  { 2.0, 3.0, 3.0, 15.0 },  { 5.0, -3, 1.0, 14.0 } };  gaussianElimination(mat);  return 0; } 

N = 3 def gaussianElimination(mat): for k in range(N): i_max = max(range(k, N), key=lambda i: abs(mat[i][k])) if mat[i_max][k] == 0: print("Singular Matrix.") print("Inconsistent System." if mat[i_max][N] else "May have infinitely many solutions.") return mat[k], mat[i_max] = mat[i_max], mat[k] for i in range(k + 1, N): f = mat[i][k] / mat[k][k] mat[i][k + 1:] = [mat[i][j] - f * mat[k][j] for j in range(k + 1, N + 1)] mat[i][k] = 0 x = [0] * N for i in range(N - 1, -1, -1): x[i] = (mat[i][N] - sum(mat[i][j] * x[j] for j in range(i + 1, N))) / mat[i][i] print("\nSolution for the system:") for xi in x: print(f"{xi:.8f}") mat = [[3.0, 2.0, -4.0, 3.0], [2.0, 3.0, 3.0, 15.0], [5.0, -3.0, 1.0, 14.0]] gaussianElimination(mat) 

using System; class GFG {  public static int N = 3;  static void gaussianElimination(double[, ] mat)  {  int singular_flag = forwardElim(mat);  if (singular_flag != -1) {  Console.WriteLine("Singular Matrix.");  if (mat[singular_flag, N] != 0)  Console.Write("Inconsistent System.");  else  Console.Write(  "May have infinitely many solutions.");  return;  }  backSub(mat);  }  static void swap_row(double[, ] mat, int i, int j)  {  for (int k = 0; k <= N; k++) {  double temp = mat[i, k];  mat[i, k] = mat[j, k];  mat[j, k] = temp;  }  }  static void print(double[, ] mat)  {  for (int i = 0; i < N; i++, Console.WriteLine())  for (int j = 0; j <= N; j++)  Console.Write(mat[i, j]);  Console.WriteLine();  }  static int forwardElim(double[, ] mat)  {  for (int k = 0; k < N; k++) {  int i_max = k;  int v_max = (int)mat[i_max, k];  for (int i = k + 1; i < N; i++) {  if (Math.Abs(mat[i, k]) > v_max) {  v_max = (int)mat[i, k];  i_max = i;  }  if (mat[k, i_max] == 0)  return k;  if (i_max != k)  swap_row(mat, k, i_max);  for (int i = k + 1; i < N; i++) {  double f = mat[i, k] / mat[k, k];  for (int j = k + 1; j <= N; j++)  mat[i, j] -= mat[k, j] * f;  mat[i, k] = 0;  }  }  }  return -1;  }  static void backSub(double[, ] mat)  {  double[] x = new double[N];  for (int i = N - 1; i >= 0; i--) {  x[i] = mat[i, N];  for (int j = i + 1; j < N; j++) {  x[i] -= mat[i, j] * x[j];  }  x[i] = x[i] / mat[i, i];  }  Console.WriteLine();  Console.WriteLine("Solution for the system:");  for (int i = 0; i < N; i++) {  Console.Write("{0:F6}", x[i]);  Console.WriteLine();  }  }  public static void Main(String[] args)  {  double[, ] mat = { { 3.0, 2.0, -4.0, 3.0 },  { 2.0, 3.0, 3.0, 15.0 },  { 5.0, -3, 1.0, 14.0 } };  gaussianElimination(mat);  } } 

let N = 3; function gaussianElimination(mat) {  let singular_flag = forwardElim(mat);  if (singular_flag != -1) {  console.log("Singular Matrix.");  if (mat[singular_flag][N])  console.log("Inconsistent System.");  else  console.log(  "May have infinitely many solutions.");  return;  }  backSub(mat); } function swap_row(mat, i, j) {  for (var k = 0; k <= N; k++) {  let temp = mat[i][k];  mat[i][k] = mat[j][k];  mat[j][k] = temp;  } } function print(mat) {  for (var i = 0; i < N; i++, console.log(""))  for (var j = 0; j <= N; j++)  process.stdout.write("" + mat[i][j]);  console.log(""); } function forwardElim(mat) {  for (var k = 0; k < N; k++) {  var i_max = k;  var v_max = mat[i_max][k];  for (var i = k + 1; i < N; i++)  if (Math.abs(mat[i][k]) > v_max)  v_max = mat[i][k], i_max = i;  if (!mat[k][i_max])  return k;  if (i_max != k)  swap_row(mat, k, i_max);  for (var i = k + 1; i < N; i++) {  let f = mat[i][k] / mat[k][k];  for (var j = k + 1; j <= N; j++)  mat[i][j] -= mat[k][j] * f;  mat[i][k] = 0;  }  }  return -1; } function backSub(mat) {  let x = new Array(N);  for (var i = N - 1; i >= 0; i--) {  x[i] = mat[i][N];  for (var j = i + 1; j < N; j++) {  x[i] -= mat[i][j] * x[j];  }  x[i] = Math.round(x[i] / mat[i][i]);  }  console.log("\nSolution for the system:");  for (var i = 0; i < N; i++)  console.log(x[i].toFixed(8)); } let mat = [  [ 3.0, 2.0, -4.0, 3.0 ], [ 2.0, 3.0, 3.0, 15.0 ],  [ 5.0, -3, 1.0, 14.0 ] ]; gaussianElimination(mat); 

Solution for the system: 3.000000 1.000000 2.000000 

Illustration: 

2. Partial Pivoting

The steps involved in Partial Pivoting Gaussian Elimination are:

1. Start with the given augmented matrix.
2. Find the row with the largest absolute value in the first column and swap that row with the first row.
3. Divide the first row by the pivot element to make it equal to 1.
4. Use the first row to eliminate the first column in all other rows by subtracting the appropriate multiple of the first row from each row.
5. Repeat steps 2-4 for the remaining columns, using the row with the largest absolute value as the pivot row for each column.
6. Back-substitute to obtain the values of the unknowns.

#include <iostream> #include <cmath> using namespace std; const int MAXN = 100; void partial_pivot(double A[MAXN][MAXN+1], int n) {  for (int i = 0; i < n; i++) {  int pivot_row = i;  for (int j = i+1; j < n; j++) {  if (abs(A[j][i]) > abs(A[pivot_row][i])) {  pivot_row = j;  }  }  if (pivot_row != i) {  for (int j = i; j <= n; j++) {  swap(A[i][j], A[pivot_row][j]);  }  }  for (int j = i+1; j < n; j++) {  double factor = A[j][i] / A[i][i];  for (int k = i; k <= n; k++) {  A[j][k] -= factor * A[i][k];  }  }  } } void back_substitute(double A[MAXN][MAXN+1], int n, double x[MAXN]) {  for (int i = n-1; i >= 0; i--) {  double sum = 0;  for (int j = i+1; j < n; j++) {  sum += A[i][j] * x[j];  }  x[i] = (A[i][n] - sum) / A[i][i];  } } int main() {  int n = 3;  double A[MAXN][MAXN+1] = {{3.0, 2.0,-4.0, 3.0},  {2.0, 3.0, 3.0, 15.0},  {5.0, -3, 1.0, 14.0}  };  double x[MAXN];    partial_pivot(A, n);  back_substitute(A, n, x);    cout << "Solution for the system:\n";  for (int i = 0; i < n; i++) {    cout << x[i] << endl;  }    return 0; } 

public class Main {    private static final int MAXN = 100;    public static void partialPivot(double[][] A, int n) {  for (int i = 0; i < n; i++) {  int pivotRow = i;  for (int j = i+1; j < n; j++) {  if (Math.abs(A[j][i]) > Math.abs(A[pivotRow][i])) {  pivotRow = j;  }  }  if (pivotRow != i) {  for (int j = i; j <= n; j++) {  double temp = A[i][j];  A[i][j] = A[pivotRow][j];  A[pivotRow][j] = temp;  }  }  for (int j = i+1; j < n; j++) {  double factor = A[j][i] / A[i][i];  for (int k = i; k <= n; k++) {  A[j][k] -= factor * A[i][k];  }  }  }  }    public static void backSubstitute(double[][] A, int n, double[] x) {  for (int i = n-1; i >= 0; i--) {  double sum = 0;  for (int j = i+1; j < n; j++) {  sum += A[i][j] * x[j];  }  x[i] = (A[i][n] - sum) / A[i][i];  }  }    public static void main(String[] args) {  int n = 3;  double[][] A = {{3.0, 2.0,-4.0, 3.0},  {2.0, 3.0, 3.0, 15.0},  {5.0, -3, 1.0, 14.0}  };  double[] x = new double[MAXN];    partialPivot(A, n);  backSubstitute(A, n, x);    System.out.println("Solution for the system:");  for (int i = 0; i < n; i++) {  System.out.println(x[i]);  }  } } 

import numpy as np MAXN = 100 def partial_pivot(A, n): for i in range(n): pivot_row = i for j in range(i + 1, n): if abs(A[j][i]) > abs(A[pivot_row][i]): pivot_row = j if pivot_row != i: A[[i, pivot_row]] = A[[pivot_row, i]] for j in range(i + 1, n): factor = A[j][i] / A[i][i] A[j] -= factor * A[i] def back_substitute(A, n): x = np.zeros(n) for i in range(n - 1, -1, -1): sum_val = sum(A[i][j] * x[j] for j in range(i + 1, n)) x[i] = (A[i][n] - sum_val) / A[i][i] return x if __name__ == "__main__": n = 3 A = np.array([[3.0, 2.0, -4.0, 3.0], [2.0, 3.0, 3.0, 15.0], [5.0, -3, 1.0, 14.0]]) partial_pivot(A, n) x = back_substitute(A, n) print("Solution for the system:") for i in range(n): print(f"{x[i]:.6f}") 

using System; class Program {  const int MAXN = 100;  static void partial_pivot(double[, ] A, int n)  {  for (int i = 0; i < n; i++) {  int pivot_row = i;  for (int j = i + 1; j < n; j++) {  if (Math.Abs(A[j, i])  > Math.Abs(A[pivot_row, i])) {  pivot_row = j;  }  }  if (pivot_row != i) {  for (int j = i; j <= n; j++) {  double temp = A[i, j];  A[i, j] = A[pivot_row, j];  A[pivot_row, j] = temp;  }  }  for (int j = i + 1; j < n; j++) {  double factor = A[j, i] / A[i, i];  for (int k = i; k <= n; k++) {  A[j, k] -= factor * A[i, k];  }  }  }  }  static void back_substitute(double[, ] A, int n,  double[] x)  {  for (int i = n - 1; i >= 0; i--) {  double sum = 0;  for (int j = i + 1; j < n; j++) {  sum += A[i, j] * x[j];  }  x[i] = (A[i, n] - sum) / A[i, i];  }  }  static void Main(string[] args)  {  int n = 3;  double[, ] A = new double[, ]  {  { 3.0, 2.0, -4.0, 3.0 }, 

const MAXN = 100; function partial_pivot(A, n) {  for (let i = 0; i < n; i++) {  let pivot_row = i;  for (let j = i + 1; j < n; j++) {  if (Math.abs(A[j][i])  > Math.abs(A[pivot_row][i])) {  pivot_row = j;  }  }  if (pivot_row != i) {  for (let j = i; j <= n; j++) {  [A[i][j], A[pivot_row][j]] =  [ A[pivot_row][j], A[i][j] ];  }  }  for (let j = i + 1; j < n; j++) {  let factor = A[j][i] / A[i][i];  for (let k = i; k <= n; k++) {  A[j][k] -= factor * A[i][k];  }  }  } } function back_substitute(A, n, x) {  for (let i = n - 1; i >= 0; i--) {  let sum = 0;  for (let j = i + 1; j < n; j++) {  sum += A[i][j] * x[j];  }  x[i] = (A[i][n] - sum) / A[i][i];  } } let n = 3; let A = [  [ 3.0, 2.0, -4.0, 3.0 ], [ 2.0, 3.0, 3.0, 15.0 ],  [ 5.0, -3, 1.0, 14.0 ] ]; let x = Array(n); partial_pivot(A, n); back_substitute(A, n, x); console.log("Solution for the system:"); for (let i = 0; i < n; i++) {  console.log(Math.round(x[i])); } 

Solution for the system: 3 1 2 

Applications of Gaussian Elimination Method

• Computing Determinants: It simplifies finding the determinant of a square matrix by reducing it to row echelon form. The determinant is calculated by multiplying the diagonal elements and adjusting for row swaps and scalar multiplications.

• Finding the Inverse of a Matrix: Gauss–Jordan elimination, a variant of Gaussian elimination, is used to find the inverse of a matrix. By augmenting the matrix with the identity matrix and applying row operations, the matrix is transformed into the inverse if it exists.

• Computing Ranks and Bases: It can be applied to any matrix to determine its rank and the basis for its column space. The row echelon form reveals information about the number of linearly independent rows and the columns that form a basis for the matrix's column space.

Gaussian Elimination Method Practice Problems

1. Solve the following system of equations using Gaussian elimination method.

x + y + z = 6
3x + y + 6z = 10
2x +4 y – z = 1

2. Solve the following linear system using the Gaussian elimination method.

x – y = -6
6x – 2y = 0

3. Using Gaussian elimination method, solve:

x – 9y + 3z = 1
x + y + z = 4
2x – y + z = 2