矩阵的高斯消去法（Gauss-Jordan方法）的Python实现

高斯消去法的改进形式为Gauss-Jordan Elimination Method，要求每一行的主元素所在列元素全部消去为0，除了主元素本身。区别如图：

目录：1 算法讲解 2 代码实现

代码目标：能解方阵、非方阵、给定精度的病态方程的通用Gauss-Jordan Method。

关键问题：

1 【最难的步骤】如何寻找pivot元素：自左向右，自上向下，寻找首个非0的元素，圈起来。保证自上向下每一行都有pivot元素，如果是0，就向下找同列不为0的一行，和当前行交换。

2 pivot所在行除以pivot值，令pivot为1

3 然后将pivot所在列全部消为0，效果如下图。

4 然后循环该过程，直到每一列都消除完毕

 

代码实现如下：

# -*- coding: utf-8 -*-
# @Author : ZhaoKe
# @Time : 2022-09-05 23:34
from typing import List


# input a augmented matrix, output its simpler form
def GaussJordanMethod(matrix: List[List], prec=3) -> List[List]:
    m, n = len(matrix), len(matrix[0])
    pi_r, pi_c = 0, 0
    # from left to right
    for j in range(n):
        # 关键在于寻找pivot
        # 主元素为0的话，要向下寻找同列非零的行，交换行，令这一行主元素不为0
        pi_c = j
        pivot = matrix[pi_r][pi_c]
        if pivot == 0:
            for k in range(pi_r+1, m):
                if matrix[k][pi_c] == 0:
                    continue
                else:
                    matrix[pi_r], matrix[k] = matrix[k], matrix[pi_r]
                    pivot = matrix[pi_r][j]
                    print("exchange ", matrix)
                    break
        if pivot == 0:
            continue
        # from above to bottom
        # 为了实现Gauss-Jordan，pivot位置的元素应该置为1
        print("row", pi_r, matrix[pi_r])
        if pivot != 1:
            for l in range(0, n):
                matrix[pi_r][l] = matrix[pi_r][l] / pivot
        print("row", pi_r, matrix[pi_r])
        # 不再拘泥于对角线下方消除，整整一列都要消除
        for i in range(m):
            # 主元素不可以消去，直接跳过该行
            if i == pi_r:
                continue
            # 当前行的该列元素为0的话，跳过即可
            if matrix[i][pi_c] == 0:
                continue
            # 初等变换
            # replace the jth equation by a combination of itself plus a multiple of the ith equation
            coef = matrix[i][pi_c]  # / matrix[j][j]  matrix[j][j]必为1（根据Gauss-Jordan法）
            print("coef:", coef)
            for k in range(n):
                matrix[i][k] = round(matrix[i][k] - coef * matrix[pi_r][k], prec)
            print("-> ", matrix[i])
        pi_r += 1
        if pi_r >= m:
            break
    # elimination end
    print(matrix)

    # 解方程，如果有唯一解的话
    # 如果无解，就没意义；如果有无数多个解，那么找出最大线性无关组即可（对每一行的非零元素计数然后判断就行）
    # # solution as follow
    # for i in range(m - 1, -1, -1):
    #     for j in range(n - 1, -1, -1):
    #         if matrix[i][j] == 0:
    #             continue
    #         else:
    #             print("x" + str(i+1), "=", matrix[i][-1])
    #             break
    return matrix


if __name__ == '__main__':
    # # 非方阵
    # input_m5 = [[1, 2, 1, 1], [2, 4, 2, 2], [3, 6, 3, 4]]  # 结果正确
    # input_m5 = [[1, 2, 1, 3, 3], [2, 4, 0, 4, 4], [1, 2, 3, 5, 5], [2, 4, 0, 4, 7]]  # 结果正确
    input_m5 = [[1, 1, 2, 2, 1, 1], [2, 2, 4, 4, 3, 1], [2, 2, 4, 4, 2, 2], [3, 5, 8, 6, 5, 3]]
    # input_m5 = [[1, 2, 3], [2, 6, 8], [2, 6, 0], [1, 2, 5], [3, 8, 6]]
    GaussJordanMethod(input_m5)
    # 方阵
    # input_m0 = [[2, 1, 1, 1], [6, 2, 1, -1], [-2, 2, 1, 7]]  # 结果正确
    input_m1 = [[0, 1, -1, 3], [-2, 4, -1, 1], [-2, 5, -4, -2]]  # 结果正确
    # input_m1 = [[2, 2, 6, 4], [2, 1, 7, 6], [-2, -6, -7, -1]]  # 结果正确
    GaussJordanMethod(input_m1)
    # 病态方程：
    # input_m4 = [[47, 28, 19], [89, 53, 36]]  # 结果正确
    # input_m4 = [[0.835, 0.667, 0.168], [0.333, 0.266, 0.067]]  # 结果正确
    # GaussJordanMethod(input_m4, prec=5)
    # print("=================")
    input_m4 = [[0.835, 0.667, 0.168], [0.333, 0.266, 0.067]]  # 结果正确
    GaussJordanMethod(input_m4, prec=6)

 

有效数字为5的时候，经过高斯消去的矩阵为：

[[1.0, 0.798802395209581, 0.20119760479041918], [0.0, -0.0, 0.0]]

显然无解

有效数字为6的时候，消去得：

[[1.0, 0.0, 1.0], [-0.0, 1.0, -1.0]]
x2 = -1.0
x1 = 1.0