Use gaussian elimination to solve the system of equations

View Discussion

Improve Article

Save Article

View Discussion

Improve Article

Save Article

The article focuses on using an algorithm for solving a system of linear equations. We will deal with the matrix of coefficients. Gaussian Elimination does not work on singular matrices (they lead to division by zero).

Row echelon form: Matrix is said to be in r.e.f. if the following conditions hold: 

1. The first non-zero element in each row, called the leading coefficient, is 1.
2. Each leading coefficient is in a column to the right of the previous row leading coefficient.
3. Rows with all zeros are below rows with at least one non-zero element.

Reduced row echelon form: Matrix is said to be in r.r.e.f. if the following conditions hold – 

1. All the conditions for r.e.f.
2. The leading coefficient in each row is the only non-zero entry in its column.

The algorithm is majorly about performing a sequence of operations on the rows of the matrix. What we would like to keep in mind while performing these operations is that we want to convert the matrix into an upper triangular matrix in row echelon form. The operations can be: 

1. Swapping two rows
2. Multiplying a row by a non-zero scalar
3. Adding to one row a multiple of another

The process:  

1. Forward elimination: reduction to row echelon form. Using it one can tell whether there are no solutions, or unique solution, or infinitely many solutions.
2. Back substitution: further reduction to reduced row echelon form.

Algorithm:  

1. Partial pivoting: Find the kth pivot by swapping rows, to move the entry with the largest absolute value to the pivot position. This imparts computational stability to the algorithm.
2. For each row below the pivot, calculate the factor f which makes the kth entry zero, and for every element in the row subtract the fth multiple of the corresponding element in the kth row.
3. Repeat above steps for each unknown. We will be left with a partial r.e.f. matrix.

Below is the implementation of above algorithm.  

C++

#include<bits/stdc++.h>

using namespace std;

#define N 3        // Number of unknowns

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;

}

Java

import java.io.*;

class GFG

{

  public static int N = 3;

  static void gaussianElimination(double mat[][])

  {

    int singular_flag = forwardElim(mat);

    if (singular_flag != -1)

    {

      System.out.println("Singular Matrix.");

      if (mat[singular_flag][N] != 0)

        System.out.print("Inconsistent System.");

      else

        System.out.print(

        "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++, System.out.println())

      for (int j = 0; j <= N; j++)

        System.out.print(mat[i][j]);

    System.out.println();

  }

  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];

    }

    System.out.println();

    System.out.println("Solution for the system:");

    for (int i = 0; i < N; i++)

    {

      System.out.format("%.6f", x[i]);

      System.out.println();

    }

  }

  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);

  }

}

Python3

N = 3

def gaussianElimination(mat):

    singular_flag = forwardElim(mat)

    if (singular_flag != -1):

        print("Singular Matrix.")

        if (mat[singular_flag][N]):

            print("Inconsistent System.")

        else:

            print("May have infinitely many solutions.")

        return

    backSub(mat)

def swap_row(mat, i, j):

    for k in range(N + 1):

        temp = mat[i][k]

        mat[i][k] = mat[j][k]

        mat[j][k] = temp

def forwardElim(mat):

    for k in range(N):

        i_max = k

        v_max = mat[i_max][k]

        for i in range(k + 1, N):

            if (abs(mat[i][k]) > v_max):

                v_max = mat[i][k]

                i_max = i

        if not mat[k][i_max]:

            return k   

        if (i_max != k):

            swap_row(mat, k, i_max)

        for i in range(k + 1, N):

            f = mat[i][k]/mat[k][k]

            for j in range(k + 1, N + 1):

                mat[i][j] -= mat[k][j]*f

            mat[i][k] = 0

    return -1

def backSub(mat):

    x = [None for _ in range(N)]   

    for i in range(N-1, -1, -1):

        x[i] = mat[i][N]

        for j in range(i + 1, N):

            x[i] -= mat[i][j]*x[j]

        x[i] = (x[i]/mat[i][i])

    print("\nSolution for the system:")

    for i in range(N):

        print("{:.8f}".format(x[i]))

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)

C#

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);

}

}

PHP

<?php

$N = 3;

function gaussianElimination($mat)

{

    global $N;

    $singular_flag = forwardElim($mat);

    if ($singular_flag != -1)

    {

        print("Singular Matrix.\n");

        if ($mat[$singular_flag][$N])

            print("Inconsistent System.");

        else

            print("May have infinitely many solutions.");

        return;

    }

    backSub($mat);

}

function swap_row(&$mat, $i, $j)

{

    global $N;

    for ($k = 0; $k <= $N; $k++)

    {

        $temp = $mat[$i][$k];

        $mat[$i][$k] = $mat[$j][$k];

        $mat[$j][$k] = $temp;

    }

}

function print1($mat)

{

    global $N;

    for ($i=0; $i<$N; $i++, print("\n"))

        for ($j=0; $j<=$N; $j++)

            print($mat[$i][$j]);

    print("\n");

}

function forwardElim(&$mat)

{

    global $N;

    for ($k=0; $k<$N; $k++)

    {

        $i_max = $k;

        $v_max = $mat[$i_max][$k];

        for ($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 ($i = $k + 1; $i < $N; $i++)

        {

            $f = $mat[$i][$k]/$mat[$k][$k];

            for ($j = $k + 1; $j <= $N; $j++)

                $mat[$i][$j] -= $mat[$k][$j] * $f;

            $mat[$i][$k] = 0;

        }

    }

    return -1;

}

function backSub(&$mat)

{

    global $N;

    $x = array_fill(0, $N, 0);

    for ($i = $N - 1; $i >= 0; $i--)

    {

        $x[$i] = $mat[$i][$N];

        for ($j = $i + 1; $j < $N; $j++)

        {

            $x[$i] -= $mat[$i][$j] * $x[$j];

        }

        $x[$i] = $x[$i] / $mat[$i][$i];

    }

    print("\nSolution for the system:\n");

    for ($i = 0; $i < $N; $i++)

        print(number_format(strval($x[$i]), 6)."\n");

}

    $mat = array(array(3.0, 2.0,-4.0, 3.0),

                        array(2.0, 3.0, 3.0, 15.0),

                        array(5.0, -3, 1.0, 14.0));

    gaussianElimination($mat);

?>

Javascript

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);

Output: 

Illustration: 
 

Time Complexity: Since for each pivot we traverse the part to its right for each row below it, O(n)*(O(n)*O(n)) = O(n3).
We can also apply Gaussian Elimination for calculating:

1. Rank of a matrix
2. Determinant of a matrix
3. Inverse of an invertible square matrix

This article is contributed by Yash Varyani. Please write comments if you find anything incorrect, or you want to share more information about the topic discussed above.