Solving a system of linear equations using Python to implement the elimination method (Gaussian elimination, column pivot, Gaussian Jordan)

1. Gaussian elimination method

Gaussian elimination method is a basic method for solving linear equations. Its main idea is to convert the original coefficient matrix into an upper triangular matrix or lower triangular matrix through matrix transformation, and then solve the linear equations through back substitution and other methods. solution.

Specifically, the Gaussian elimination method has the following steps:

Construct an augmented matrix: combine the coefficient matrix and the constant vector into an augmented matrix.
Forward elimination: Convert the augmented matrix into an upper triangular matrix or a lower triangular matrix through a series of elementary row transformations. The main operation in the process is addition and subtraction, that is, using the multiple of row ii to add to row jj (i< ji
Back-substitution solution: Starting from the last row, find the values of the unknowns one by one.

It should be noted that if all elements of a certain row are 00 during the forward elimination process, the elimination operation cannot be performed. At this time, a row exchange operation is required so that the row is not 00.

The advantage of Gaussian elimination method is that it is simple to understand, easy to implement, and can handle general linear equations. However, there are also some problems. Due to the impact of rounding errors and zero elements during the operation, the calculation results may be unstable for linear equations with ill-conditioned matrices or large coefficient matrices.

$\left\{ \begin{array}{l} 9{x_1} - {x_2} - {x_3} = 7\ - {x_1} + 8{x_2} = 7 \ - {x_1} + 9{x_3} = 8 \end{array} \right.$

"""
@Time: 2023/10/19 0019 17:11
@Auth:yeqc
"""
import numpy as np


def gaussian_elimination(A, b):
    """
    Gaussian elimination method to solve linear equations Ax = b

    parameter:
    A: coefficient matrix, shape (n, n)
    b: constant vector on the right side with shape (n,)

    return:
    x: solution vector with shape (n,)
    """
    n = len(b)

    # Combine the coefficient matrix and the constant on the right
    Ab = np.column_stack((A.astype(float), b.astype(float)))

    # Forward elimination
    for i in range(n - 1):
        for j in range(i + 1, n):
            factor = Ab[j, i] / Ab[i, i]
            Ab[j, i:] -= factor * Ab[i, i:]

    # Back to the solution
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, n] - np.dot(Ab[i, i:n], x[i:n])) / Ab[i, i]

    return x


# Example usage
A = np.array([[9,-1,-1],[-1,8,0],[-1,0,9]])

b = np.array([7,7,8])

x = gaussian_elimination(A, b)
print("Solution x =", x)

2. Column pivot elimination method

Column Principal Pivot Elimination Method, also known as column principal elimination method, is a method for solving linear equations. It is an improved version of the Gauss-Jordan elimination method and is used to transform augmented matrices into their row simplest form.

The basic idea of the column pivot elimination method is to select the main element of each column (also called column pivot), place the main element at the top of the corresponding column through a row exchange operation, and then use a series of row operations to place the main element at the top of the corresponding column. Other elements are eliminated.

Specific steps are as follows:

Construct an augmented matrix: combine the coefficient matrix and the constant vector into an augmented matrix.
Select column pivot: For each column, select the element with the largest absolute value in the column as the pivot element in the remaining rows, and exchange its row with the currently processed row.
Forward elimination: Convert the augmented matrix into the simplest row form through a series of elementary row transformations, even if other elements in the column of the main element become 0.
Back-substitution solution: Starting from the last row, find the values of the unknowns one by one.

Using the column pivot elimination method can avoid the numerical instability problems that may occur in the Gaussian elimination method and obtain a more accurate solution. However, compared with the Gaussian elimination method, the column pivot elimination method requires a larger amount of calculations because the largest principal element needs to be selected for row exchange in each step, which increases the complexity of the calculation. But for some special linear equations, the column-pivot elimination method may be more effective.

"""
@Time: 2023/10/19 0019 17:11
@Auth:yeqc
"""
import numpy as np


def gaussian_elimination(A, b):
    n = len(A)

    #Convert the data type of the array to float64
    A = A.astype(np.float64)
    b = b.astype(np.float64)

    # Gaussian elimination
    for i in range(n - 1):
        max_idx = i

        # Select column pivot
        for j in range(i + 1, n):
            if abs(A[j][i]) > abs(A[max_idx][i]):
                max_idx = j

        # swap lines
        A[[i, max_idx]] = A[[max_idx, i]]
        b[[i, max_idx]] = b[[max_idx, i]]

        for j in range(i + 1, n):
            # Calculate multiples
            multiplier = A[j][i] / A[i][i]

            # Update matrix
            A[j][i:] -= multiplier * A[i][i:]
            b[j] -= multiplier * b[i]

    # Back to the solution
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - np.dot(A[i][i + 1:], x[i + 1:])) / A[i][i]

    return x


# Example
A = np.array([[9, -1, -1], [-1, 8, 0], [-1, 0, 9]])

b = np.array([7, 7, 8])

x = gaussian_elimination(A, b)
print("Solution:", x)

3. Gauss-Jordan elimination method (Gauss-Jordan)

Gauss-Jordan Elimination Method is a basic method for solving linear equations. It transforms the augmented matrix into its simplest row form through a series of elementary row transformations to obtain solutions to the linear equations.

Specific steps are as follows:

Construct an augmented matrix: combine the coefficient matrix and the constant vector into an augmented matrix.
Forward elimination: Convert the augmented matrix into an upper triangular form through a series of elementary row transformations, that is, a form in which the lower left corner of the coefficient matrix is all 0.
Reverse elimination: The augmented matrix is further transformed into a diagonal form through a series of elementary row transformations, that is, a form in which all elements of the coefficient matrix except the diagonal are 0. The main operation in the process is to use addition, subtraction and elimination again to make the coefficient of each unknown quantity 1, and at the same time eliminate the coefficients of other rows to 0.
Solving the system of equations: The solution of the system of linear equations can be easily found based on the augmented matrix in diagonal form.

It should be noted that in the process of converting to upper triangular form, there may be cases where the principal element is 0 or a smaller value, causing numerical instability in the algorithm. Therefore, in the Gauss-Jordan elimination method, it is necessary to strengthen the control of numerical accuracy to avoid abnormal results.

"""
@Time: 2023/10/19 0019 17:11
@Auth:yeqc
"""
import numpy as np


def gaussian_jordan_elimination(A, b):
    n = len(A)

    #Construct augmented matrix
    augmented_matrix = np.column_stack((A.astype(float), b.astype(float)))

    # Gaussian Jordan Elimination
    for i in range(n):
        pivot = augmented_matrix[i][i]

        # Normalize the pivot to 1
        augmented_matrix[i] /= pivot

        # Update other rows
        for j in range(n):
            if j != i:
                multiplier = augmented_matrix[j][i]
                augmented_matrix[j] -= multiplier * augmented_matrix[i]

    # Extract solution
    x = augmented_matrix[:, n]

    return x


# Example
A = np.array([[9, -1, -1], [-1, 8, 0], [-1, 0, 9]])

b = np.array([7, 7, 8])

x = gaussian_jordan_elimination(A, b)
print("Solution:", x)