python实现高斯消元法求线性方程组的解

高斯消元法简介

  数学上，高斯消元法（或译：高斯消去法），是线性代数规划中的一个算法，可用来为线性方程组求解。但其算法十分复杂，不常用于加减消元法，求出矩阵的秩，以及求出可逆方阵的逆矩阵。不过，如果有过百万条等式时，这个算法会十分省时。一些极大的方程组通常会用迭代法以及花式消元来解决。当用于一个矩阵时，高斯消元法会产生出一个“行梯阵式”。高斯消元法可以用在电脑中来解决数千条等式及未知数。亦有一些方法特地用来解决一些有特别排列的系数的方程组。        ——转自百度百科

内容

  消元法是将方程组中的一方程的未知数用含有另一未知数的代数式表示，并将其代入到另一方程中，这就消去了一未知数，得到一解；或将方程组中的一方程倍乘某个常数加到另外一方程中去，也可达到消去一未知数的目的。消元法主要用于二元一次方程组的求解。        ——转自百度百科

核心

 1)两方程互换，解不变；
 2)一方程乘以非零数k，解不变；
 3)一方程乘以数k加上另一方程，解不变        ——转自百度百科

算法思路

  思路就是根据高斯消元法的特性，先将矩阵化成行阶梯型矩阵，由行阶梯型矩阵得出秩再判断线性方程组解的个数
    如果R(A)<=R(Ab) 则方程无解
    R(A)=R(Ab)=N时方程有唯一解 （齐次方程称无非0解）
    R(A)=R(Ab)<N时方程有无穷个解 （齐次方程称有非0解）
得出解的个数后，如果存在唯一解，再将之前的矩阵化为行最简形矩阵，得出线性方程组的解

算法

import numpy as np

def swag(a, b):     #交换矩阵某两行函数
    for i in range(0, len(a)):
        t = a[i]
        a[i] = b[i]
        b[i] = t

def print_matrix( info, matrix ):       #打印矩阵
    print (info)
    for i in range( 0, matrix.shape[0]):
        print('[', end='')
        for j in range( 0, matrix.shape[1]):
            if(j == matrix.shape[1] - 1):   # 向量 b
                print( '|', end=''),                           # 将系数矩阵和向量b用“|”隔开
            print("%5.2f" %matrix[i][j], end=' ')
            if j == matrix.shape[1] - 1:                  # 输出矩阵元素m[i][j]
                print(']', end=' ')
                print('\n')

def check(matrix, i, row, col):         #判定函数，用于检测当前处理的矩阵的行的值是否全为0，若全为零，则从矩阵最底行开始遍历找出不全为零的行与此行交换，再对交换后的行进行处理
    if 0.00 in set(matrix[i]) and len(set(matrix[i])) == 1:     #判定当前行是否全为零，将当前行集合化看其是否只有一个0元素
        for j in range(row - 1, i ,-1):                         #从矩阵最后一行开始寻找不全为零的一行
            try:
                if not(0.00 in set(matrix[j]) and len(set(matrix[j])) == 1):    #判定这行是否全为零
                    swag(matrix[i], matrix[j])                                  #交换
                    select(matrix, i, col)                                      #交换后再进行阶梯化处理
                    break
            except:                     #如果遍历发现没有不为零的行了，则直接返回
                return

def select(matrix, i, col):             #矩阵的阶梯化处理函数
    if 0.00 in set(matrix[i]) and len(set(matrix[i])) == 1:     #主要针对矩阵最后一行，judge函数无法对最后一行进行非零判定
        return
    for k in range(0, i):                                       #i是当前处理的矩阵的行数，这个循环是对矩阵的第i行进行阶梯化处理
        temp = matrix[i][k] / matrix[k][k]                      #根据主元matrix[k][k]确定让matrix[i][k]归零的除数
        if temp == 0:                                           #主元为0不执行操作
            continue
        for j in range(0, col):                                 #对这一行的所有数进行操作，其中和主元同一列的那个数会归零
            matrix[i][j] = matrix[i][j] - matrix[k][j] * temp

def solve(matrix):                      #算法的主体部分
    row = matrix.shape[0]                       #得到矩阵的行数
    col = matrix.shape[1]                       #得到矩阵的列数
    for i in range(0, row):                     #检查矩阵某一行的主元是否为0，为零则遍历寻找不为零的行与其交换
        if matrix[i][i] == 0:                   #易得matrix[i][i]就为第i行主元
            for j in range(i + 1, row):
                if matrix[j][i] != 0:
                    swag(matrix[i], matrix[j])
                    break
        select(matrix, i, col)                  #进行阶梯化操作
        check(matrix, i, row, col)              #判定阶梯化后会不会使此行的值全化为零，若全化为零则将其与其他不为零的行交换，再进行阶梯化处理

def to_one(matrix):                     #将矩阵化为行最简形矩阵
    row = matrix.shape[0]
    col = matrix.shape[1]
    for i in range(0, row):                         #将矩阵每一行主元全部归一
        temp = matrix[i][i]
        for j in range(i, col):
            matrix[i][j] = matrix[i][j] / temp
    for i in range(0, row - 1):                     #对矩阵每一行非主元进行归零处理
        for j in range(i + 1, col - 1):             #从主元后开始遍历对每一个数进行归零处理
            temp = matrix[i][j]                     #之前将矩阵的主元都化为1，则此处令其归零的除数就是其本身 （x / 1 = x）
            for k in range(j, col):                 # matrix[i][k]是第i行的第k个值，对此行的第j个数字进行归零处理，此行其他的数字也要做相同操作
                matrix[i][k] = matrix[i][k] - matrix[j][k] * temp   #matrix[j][k]是matrix[i][k]对应的主元行上对应的列值

def judge(matrix):                      #由阶梯化后的矩阵判定线性方程组解的个数
    row = matrix.shape[0]
    col = matrix.shape[1]
    vanumlist = []
    for i in range(0, col):             #因为由solve函数确定的是一个行阶梯型矩阵，则仅需判定最后一行非零数字的个数。若无或者有多个，则证明多解；有一个则无解；有两个则证明有唯一解
        if matrix[row - 1][i] != 0:
            vanumlist.append(matrix[row - 1][i])    #每遍历到一个非零值就将其写入vanumlist这个list中，最后计算vanumlist的长度即可
    if len(vanumlist) == 1:
        print( '方程组无解')
    elif len(vanumlist) == 2:
        to_one(matrix)                  #若有唯一解，则将矩阵进一步化成行最简形矩阵求出解
        print_matrix('将矩阵最简阶梯化为：', matrix)
        for i in range(0, row):
            print("x%d = %4.2f" %(i,matrix[i][col - 1]), end="  ")
    else:
        print( '方程组有多解')

matrix = np.array([ [ 1, 1, 1, 1, 1,  1],
                    [ 3, 2, 1, 1,-3, 25],
                    [ 0, 1, 2, 2, 6,-22],
                    [ 5, 4, 3, 3,-1, 27],
                    [ 2,-1, 0, 3, 2,  2]
                    ],dtype=float)
# matrix = np.array([ [1,-2,3,-1,-1,2],
#                     [1,1,-1,1,-2,1],
#                     [2,-1,1,0,-2,1],
#                     [2,2,-5,2,-1,5],
#                     ],dtype=float)
# matrix = np.array([ [2,3,11,5,2],
#                     [1,1,5,2,1],
#                     [2,1,3,2,-3],
#                     [1,1,3,4,-3],
#                     ],dtype=float)
print_matrix('初始矩阵为：', matrix)
solve(matrix)
print_matrix('矩阵阶梯化为：', matrix)
judge(matrix)

代码执行结果