Gaussian Elimination

Linear Algebra

Gaussian elimination is a method for solving systems of linear equations by transforming the system into an upper triangular form using row operations. The process consists of forward elimination to create an upper triangular matrix, followed by back-substitution to find the solution.

To Pivot or Not to Pivot?

Gaussian elimination can sometimes be performed without pivoting, but in numerical implementations, partial pivoting is typically required. This is because roundoff errors can accumulate and cause the method to fail.

An Example

Consider the system of equations:

2 x + y - z - 3 x - y + 2 z - 2 x + y + 2 z = 8, = - 11, = - 3.

We represent this system as an augmented matrix:

2 - 3 - 2 1 - 1 1 - 1 22 8 - 11 - 3

Step 1: Normalize the First Pivot

Make the first pivot (top-left element) a 1 by dividing row 1 by 2:

1 - 3 - 2 0.5 - 1 1 - 0.5 22 4 - 11 - 3

Step 2: Eliminate Below the First Pivot

Use row 1 to eliminate the first column entries below the pivot:

Replace row 2 with $R_{2} + 3 R_{1}$ .
Replace row 3 with $R_{3} + 2 R_{1}$ .

100 0.5 0.5 2 - 0.5 0.5 1 415

Step 3: Normalize the Second Pivot

Make the second pivot (row 2, column 2) a 1 by dividing row 2 by 0.5:

100 0.5 12 - 0.5 11 425

Step 4: Eliminate Below the Second Pivot

Use row 2 to eliminate the second column entry in row 3:

Replace row 3 with $R_{3} - 2 R_{2}$ .

100 0.5 10 - 0.5 1 - 1 421

Step 5: Normalize the Last Pivot

Make the last pivot (row 3, column 3) a 1 by multiplying row 3 by -1:

100 0.5 10 - 0.5 11 42 - 1

Step 6: Back-Substitution

Now that the system is in upper triangular form, solve for the variables:

From the last row: $z = - 1$ .
From the second row: $y + z = 2 \Rightarrow y = 3$ .
From the first row: $x + 0.5 y - 0.5 z = 4$ .

Substituting known values:

x + 1.5 + 0.5 = 4 \Rightarrow x = 2

Thus, the solution to the system is:

x = 2, y = 3, z = - 1

Algorithm Implementation

Here's a simple C implementation of the Gaussian Elimination with pivoting.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAX 10

// Function to swap two rows in the matrix
void swapRows(double matrix[MAX][MAX + 1], int row1, int row2, int n) {
    for (int i = 0; i <= n; i++) {
        double temp = matrix[row1][i];
        matrix[row1][i] = matrix[row2][i];
        matrix[row2][i] = temp;
    }
}

// Gaussian elimination with partial pivoting
void gaussianElimination(double matrix[MAX][MAX + 1], int m, int n) {
    int h = 0;  // Pivot row index
    int k = 0;  // Pivot column index

    while (h < m && k < n) {
        // Find the k-th pivot: maximum absolute value in column k
        int i_max = h;
        for (int i = h + 1; i < m; i++) {
            if (fabs(matrix[i][k]) > fabs(matrix[i_max][k])) {
                i_max = i;
            }
        }

        // If pivot is zero, skip the column
        if (matrix[i_max][k] == 0) {
            k++;
        } else {
            // Swap rows
            swapRows(matrix, h, i_max, n);

            // Row reduction: for all rows below the pivot row
            for (int i = h + 1; i < m; i++) {
                double f = matrix[i][k] / matrix[h][k];
                matrix[i][k] = 0;  // Fill with zeros the lower part of pivot column

                // Update remaining elements in the row
                for (int j = k + 1; j < n; j++) {
                    matrix[i][j] = matrix[i][j] - matrix[h][j] * f;
                }
            }

            // Move to the next pivot row and column
            h++;
            k++;
        }
    }
}

// Back substitution to find the solution
void backSubstitution(double matrix[MAX][MAX + 1], double solution[MAX], int n) {
    for (int i = n - 1; i >= 0; i--) {
        solution[i] = matrix[i][n];
        for (int j = i + 1; j < n; j++) {
            solution[i] -= matrix[i][j] * solution[j];
        }
        solution[i] /= matrix[i][i];
    }
}

int main() {
    int m, n;

    printf("Enter the number of equations (m): ");
    scanf("%d", &m);
    printf("Enter the number of variables (n): ");
    scanf("%d", &n);

    double matrix[MAX][MAX + 1];
    double solution[MAX];

    printf("Enter the augmented matrix (coefficients and constants):\n");
    for (int i = 0; i < m; i++) {
        for (int j = 0; j <= n; j++) {
            scanf("%lf", &matrix[i][j]);
        }
    }

    // Perform Gaussian elimination
    gaussianElimination(matrix, m, n);

    // Perform back substitution
    backSubstitution(matrix, solution, n);

    // Output the solution
    printf("The solution is:\n");
    for (int i = 0; i < n; i++) {
        printf("x%d = %.2lf\n", i + 1, solution[i]);
    }

    return 0;
}