python——高斯消元法与高斯列主消元

        刚做完8个大作业，结束了上课生涯，十分激动，今后的学习生涯可能会困难重重，遇到困难尝尝使我怀疑人生，但只要有足够的毅力去战胜困难，像拥有科比渴望赢球的欲望，这都不是什么大问题，希望大家共勉。朋友之前说过，我的畏难情绪太过严重，经过最近的磨练，我感觉找到了一点学习的方法与乐趣，我会继续保持并跟进。

        高斯消元法与高斯列主消元法的课本知识小弟不再赘述，《数值计算》课本上可以找到相应内容，或者直接百度。高斯消元法的前提是系数矩阵为非奇异矩阵，我们采用是否满秩来等效判断，若满秩，则为非奇异矩阵。在高斯消元过程中，如果遇到主元为0或者小量时，则高斯消元不再适用，应采用高斯列主消元，其核心思想是每次消元前将绝对值最大的行与主元行进行交换。代码共分为四部分：(1)系数矩阵非奇异判断；(2)高斯消元；(3)高斯列主消元；(4)代码测试与结果。

       部分内容参考了同窗Zhou·X·G的思想(如方程回代部分，绞尽脑汁，可能是因为初次写逻辑性如此强的代码，我相信今后我的逻辑思维会更加敏捷)，他一直是我学习的榜样，在他身上我看到了学者应该的态度与热情，在此对其表示由衷的感谢。最后欢迎各位前辈对小弟的代码进行批评指正，初来乍到，还有很多需要向各位前辈们学习。

#-*- coding:utf-8 -*-
#May 24th 2021, on Mon
#Luo·Y·T

###########################
#高斯消元与高斯列主消元解方程
#不使用任何第三方库
###########################
#判断系数矩阵是否为非奇异矩阵
def IsItNonSingular(A):
    row = len(A) #系数矩阵行数
    col  = len(A[0]) #系数矩阵列数
    #控制第i步，A矩阵为n阶方阵
    for i in range(row-1): 
        #若主元元素为0，则进行行变换
        if A[i][i] == 0:
            #从后往前找出主元为0的列中，首个不为零的行
            #对主元为0的行与此行进行行变换
            Str = [] 
            for h in range(row-1, i, -1): #添加该列所有元素
                Str.append(A[h][i])
            Num = col - Str.index(next(filter(lambda x: x!=0, Str))) 
            A[Num-1], A[i] = A[i], A[Num-1] #做行变换  
        #消元运算
        for j in range(i+1, row, 1): 
            coeff = A[j][i] / A[i][i]
            for k in range(i, col, 1): 
                A[j][k] = A[j][k] - coeff * A[i][k] #系数矩阵消元
    #判断该方程组系数矩阵是否为非奇异矩阵
    for i in range(row):
        if A[i][i] == 0:
            print("Coefficient matrix is not a nonsingular matrix.")
            return 'N'
    print("Coefficient matrix is a nonsingular matrix.")

def GuassElimination(A, b):
    row = len(A) #系数矩阵行数
    col  = len(A[0]) #系数矩阵列数
    ε = 1E-5 #定义一个小量
    print("系数矩阵：", A)
    print("常数列：", b)
    #控制第i步，高斯消元需要n-1步，A矩阵为n阶方阵
    for i in range(row-1):
        if abs(A[i][i]) <= ε: #若主元为一个小量，则采用列主消元
            return None
        #消元运算
        for j in range(i+1, row, 1): 
            coeff = A[j][i] / A[i][i]
            for k in range(i, col, 1): 
                A[j][k] = A[j][k] - coeff*A[i][k] #系数矩阵消元
            b[j] = b[j] - coeff * b[i] #对应常数列消元

    #回代过程
    x=[0] * col #初始化元组，用于后面存放解
    x[col-1] = b[col-1] / A[col-1][col-1] #第n个解
    for i in range(row-2, -1, -1):
        for j in range(col-1, i, -1):
            b[i] = b[i] - A[i][j] * x[j]
        x[i] = b[i] / A[i][i]
    print("方程组的解为：")
    for i in range(col):
        print("x{}:".format(i+1), '%.8f' %x[i])
    print("\n")
    return x

def ColumnPrincipalElimination(A, b):
    row = len(A) #系数矩阵行数
    col  = len(A[0]) #系数矩阵列数
    print("系数矩阵：", A)
    print("常数列：", b)
    #控制第i步，消元需要n-1步，A矩阵为n阶方阵
    for i in range(row-1): 
        #每步运算前找出列中绝对值最大元素
        #作行变换，让绝对值最大的元素行作主元
        Str = []
        for j in range(i, row, 1):
            Str.append(A[j][i])
        Num = Str.index(max(Str)) + i
        A[Num], A[i] = A[i], A[Num] #行变换
        b[Num], b[i] = b[i], b[Num]
        
       #消元运算
        for j in range(i+1, row, 1): 
            coeff = A[j][i] / A[i][i]
            for k in range(i, col, 1): 
                A[j][k] = A[j][k] - coeff * A[i][k] #系数矩阵消元
            b[j] = b[j] - coeff * b[i] #对应常数列消元
    
    #回代过程
    x = [0] * col #初始化元组，用于后面存放解
    x[col-1] = b[col-1] / A[col-1][col-1] #第n个解
    for i in range(row-2, -1, -1):
        for j in range(col-1, i, -1):
            b[i] = b[i] - A[i][j] * x[j]
        x[i] = b[i] / A[i][i]
    print("方程组的解为：")
    for i in range(col):
        print("x{}:".format(i+1), '%.8f' %x[i])
    #print("\n")
    return x

if __name__ == '__main__':
    b1 = [1,1,1]
    b2 = [1,4,1]
    b3 = [1,12,13]
    test1 = [[1,1,1], [1,3,4], [4,5,6]]
    test2 = [[0,2,3], [4,4,4], [4,7,9]]
    test3 = [[1,1,3], [7,6,4], [8,9,9]]
    #python为动态语言，
    #在判断非奇异矩阵过程中会改动初值
    #故重新赋值
    IsItNonSingular(test1)
    print("test1高斯消元测试：")
    test1 = [[1,1,1], [1,3,4], [4,5,6]]
    GuassElimination(test1, b1)
    
    IsItNonSingular(test3)
    print("test2高斯消元测试：")
    test3 = [[1,1,3], [7,6,4], [8,9,9]]
    GuassElimination(test3, b3)

    IsItNonSingular(test2)
    print("test3高斯列主消元测试：")
    test2 = [[0,2,3], [4,4,4], [4,7,9]]
    ColumnPrincipalElimination(test2, b2)