高斯消去法的改进形式为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
Original: https://www.cnblogs.com/zhaoke271828/p/16660295.html
Author: Viktor_Cullen
Title: 矩阵的高斯消去法(Gauss-Jordan方法)的Python实现
原创文章受到原创版权保护。转载请注明出处:https://www.johngo689.com/566017/
转载文章受原作者版权保护。转载请注明原作者出处!