矩阵的高斯消去法(Gauss-Jordan方法)的Python实现

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

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

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

矩阵的高斯消去法(Gauss-Jordan方法)的Python实现

关键问题:

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

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

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

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

矩阵的高斯消去法(Gauss-Jordan方法)的Python实现

代码实现如下:

-*- 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

Original: https://www.cnblogs.com/zhaoke271828/p/16660295.html
Author: Viktor_Cullen
Title: 矩阵的高斯消去法(Gauss-Jordan方法)的Python实现

原创文章受到原创版权保护。转载请注明出处:https://www.johngo689.com/566017/

转载文章受原作者版权保护。转载请注明原作者出处!

(0)

大家都在看

亲爱的 Coder【最近整理,可免费获取】👉 最新必读书单  | 👏 面试题下载  | 🌎 免费的AI知识星球