Gaussian Elimination is a systematic method for solving systems of linear equations by transforming
the augmented matrix into row echelon form. The method is named after Carl Friedrich Gauss, although
it was known to Chinese mathematicians as early as 179 AD. It consists of a sequence of operations
performed on the corresponding matrix of coefficients.
How It Works
Forward Elimination Phase:
Start with the augmented matrix [A|b]
Convert matrix to upper triangular form
For each pivot element, eliminate all elements below it
Use row operations to make elements below pivot zero
Back Substitution Phase:
Start from the last equation
Solve for each variable
Substitute known values into previous equations
Continue until all variables are found
Pivoting Strategy:
Select largest possible pivot for numerical stability
Swap rows if necessary (partial pivoting)
Check for singular or ill-conditioned matrices
Solution Verification:
Check if solution exists
Verify solution consistency
Handle special cases (no solution, infinite solutions)
History
The method of systematically eliminating variables from equations dates back to ancient Chinese
mathematics, appearing in "The Nine Chapters on the Mathematical Art" (179 AD). However, its modern
form was developed by Carl Friedrich Gauss in the early 19th century. Gauss's work on celestial
mechanics required solving large systems of equations, leading to his refinement of the elimination
method. The technique gained prominence in the computer age due to its systematic nature and
adaptability to digital computation.
import numpy as np
def gaussian_elimination(A, b):
n = len(A)
# Augment the matrix A with the vector b
Ab = np.concatenate((A, b.reshape(n, 1)), axis=1)
# Forward elimination
for i in range(n):
# Find pivot
pivot = abs(Ab[i][i])
pivot_row = i
for k in range(i + 1, n):
if abs(Ab[k][i]) > pivot:
pivot = abs(Ab[k][i])
pivot_row = k
# Swap maximum row with current row
if pivot_row != i:
Ab[i], Ab[pivot_row] = Ab[pivot_row], Ab[i]
# Make all rows below this one 0 in current column
for k in range(i + 1, n):
factor = Ab[k][i] / Ab[i][i]
Ab[k] = Ab[k] - factor * Ab[i]
# Back substitution
x = np.zeros(n)
for i in range(n-1, -1, -1):
x[i] = Ab[i][n]
for j in range(i+1, n):
x[i] = x[i] - Ab[i][j] * x[j]
x[i] = x[i] / Ab[i][i]
return x
# Example usage
if __name__ == "__main__":
# Define system of equations
A = np.array([[2, 1, -1],
[-3, -1, 2],
[-2, 1, 2]])
b = np.array([8, -11, -3])
solution = gaussian_elimination(A, b)
print("Solution:", solution)
#include <vector>
#include <cmath>
using namespace std;
class GaussianElimination {
private:
vector<vector<double>> Ab;
int n;
void forwardElimination() {
for (int i = 0; i < n; i++) {
// Find pivot
int pivot_row = i;
double pivot_value = abs(Ab[i][i]);
for (int k = i + 1; k < n; k++) {
if (abs(Ab[k][i]) > pivot_value) {
pivot_value = abs(Ab[k][i]);
pivot_row = k;
}
}
// Swap rows if needed
if (pivot_row != i) {
swap(Ab[i], Ab[pivot_row]);
}
// Eliminate column
for (int k = i + 1; k < n; k++) {
double factor = Ab[k][i] / Ab[i][i];
for (int j = i; j <= n; j++) {
Ab[k][j] -= factor * Ab[i][j];
}
}
}
}
vector<double> backSubstitution() {
vector<double> x(n);
for (int i = n - 1; i >= 0; i--) {
x[i] = Ab[i][n];
for (int j = i + 1; j < n; j++) {
x[i] -= Ab[i][j] * x[j];
}
x[i] /= Ab[i][i];
}
return x;
}
public:
vector<double> solve(vector<vector<double>> A, vector<double> b) {
n = A.size();
Ab = vector<vector<double>>(n, vector<double>(n + 1));
// Create augmented matrix
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
Ab[i][j] = A[i][j];
}
Ab[i][n] = b[i];
}
forwardElimination();
return backSubstitution();
}
};
int main() {
vector<vector<double>> A = {
{2, 1, -1},
{-3, -1, 2},
{-2, 1, 2}
};
vector<double> b = {8, -11, -3};
GaussianElimination solver;
vector<double> solution = solver.solve(A, b);
cout << "Solution:" << endl;
for (double x : solution) {
cout << x << " ";
}
cout << endl;
return 0;
}
public class GaussianElimination {
private double[,] Ab;
private int n;
private void ForwardElimination() {
for (int i = 0; i < n; i++) {
// Find pivot
int pivotRow = i;
double pivotValue = Math.Abs(Ab[i,i]);
for (int k = i + 1; k < n; k++) {
if (Math.Abs(Ab[k,i]) > pivotValue) {
pivotValue = Math.Abs(Ab[k,i]);
pivotRow = k;
}
}
// Swap rows if needed
if (pivotRow != i) {
for (int j = i; j <= n; j++) {
double temp = Ab[i,j];
Ab[i,j] = Ab[pivotRow,j];
Ab[pivotRow,j] = temp;
}
}
// Eliminate column
for (int k = i + 1; k < n; k++) {
double factor = Ab[k,i] / Ab[i,i];
for (int j = i; j <= n; j++) {
Ab[k,j] -= factor * Ab[i,j];
}
}
}
}
private double[] BackSubstitution() {
double[] x = new double[n];
for (int i = n - 1; i >= 0; i--) {
x[i] = Ab[i,n];
for (int j = i + 1; j < n; j++) {
x[i] -= Ab[i,j] * x[j];
}
x[i] /= Ab[i,i];
}
return x;
}
public double[] Solve(double[,] A, double[] b) {
n = b.Length;
Ab = new double[n,n+1];
// Create augmented matrix
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
Ab[i,j] = A[i,j];
}
Ab[i,n] = b[i];
}
ForwardElimination();
return BackSubstitution();
}
public static void Main() {
double[,] A = new double[,] {
{2, 1, -1},
{-3, -1, 2},
{-2, 1, 2}
};
double[] b = new double[] {8, -11, -3};
var solver = new GaussianElimination();
double[] solution = solver.Solve(A, b);
Console.WriteLine("Solution:");
foreach (double x in solution) {
Console.Write($"{x:F6} ");
}
Console.WriteLine();
}
}
Complexity Analysis
Time Complexity
The time complexity is O(n³), where n is the size of the matrix:
