高斯消元法求解线性方程组（附python代码）

输入：a是m×n的系数矩阵，b是m×1的（列）向量。

输出：方程组的通解。

---

用高斯消元法（行化简法）解线性方程组步骤

• 1.构造方程组的增广矩阵
• 2.从最左边列往右，使用行化简算法把增广矩阵化为阶梯形，确定矩阵是否有解：
  
  若最后一列为主元列（最后一行非零行形如 [0 0 0 5]），无解，返回无解。
• 3.继续行化简，把主元上面的所有的元素都化为0，把主元位置变成1.
• 4.把每个主元列对应的变量表示成非主元变量的线性组合.

import numpy as np

def arguemented_mat(a, b):
    return np.c_[a,b]


# 行交换
def swap_row(a,i,j):
    m,n = a.shape
    if i >= m or j >= m:
        print 'error: out of index ...'
    else:
        a[i] = a[i] ^ a[j]
        a[j] = a[i] ^ a[j]
        a[i] = a[i] ^ a[j]
    return a


# 变成阶梯矩阵
def trape_mat(sigma):
    m,n = sigma.shape
    #保存主元所在的列数，一般来说，每行都有一个主元，除非某行全零
    main_factor = []
    main_col = 0
    while main_col < n and len(main_factor) < m:
        # 当行数多于列数的时候，出现所有的列已经处理完，结束
        if main_col == n:
            break
        # 逐列找主元,若该列全零（从第i行往下），则没有主元
        # 当前查找小矩阵首行所在原矩阵的行号
        first_row = len(main_factor)
        while main_col < n:
            new_col = sigma[first_row:, main_col]
            not_zeros = np.where(new_col > 0)[0]
            # 若该列没有非零值，则该列非主元列
            if len(not_zeros) == 0:
                main_col += 1
                break

            # 否则通过行变换找到主元
            else:
                main_factor.append(main_col)
                index = not_zeros[0]
                # 若首个元素不是主元,需要行变换
                if index != 0:
                    sigma = swap_row(sigma, first_row, first_row + index)
                # 把该主元下面的元素全部变成0
                if first_row < m - 1:
                    for k in xrange(first_row+1, m):
                        times = float(sigma[k, main_col]) / sigma[first_row, main_col]
                        sigma[k] = sigma[k] - times * sigma[first_row]
                # 处理完当前主元列之后继续
                main_col += 1
                break
    return sigma, main_factor


# 回代求解
def back_solve(sigma, main_factor):
    # 判断是否有解
    if len(main_factor) == 0:
        print 'wrong main_factor ...'
        return None
    m,n = sigma.shape
    if main_factor[-1] == n-1:
        print 'no answer ...'
        return None
    # 把所有的主元元素上方的元素变成0
    for i in range(len(main_factor)-1,-1,-1):
        factor = sigma[i, main_factor[i]]
        sigma[i] = sigma[i] / float(factor)
        for j in xrange(i):
            times = sigma[j, main_factor[i]]
            sigma[j] = sigma[j] - float(times)*sigma[i]
    # 先看看结果对不对
    return sigma

# 结果打印
def print_result(sigma, main_factor):
    if sigma is None:
        print 'no answer ...'
        return 
    m,n = sigma.shape
    result = [''] * (n-1)
    main_factor = list(main_factor)
    for i in range(n-1):
        # 如果不是主元列，则为自由变量
        if i not in main_factor:
            result[i] = 'x_' + str(i+1) + '(free var)'
        # 否则是主元变量，从对应的行，将主元变量表示成非主元变量的线性组合
        else:
            # row_of_maini表示该主元所在的行
            row_of_maini = main_factor.index(i)
            result[i] = str(sigma[row_of_maini, -1])
            for j in range(i+1, n-1):
                ratio = sigma[row_of_maini, j]
                if ratio > 0:
                    result[i] = result[i] + '-' + str(ratio) + 'x_' + str(j+1)
                if ratio < 0:
                    result[i] = result[i] + '+' + str(-ratio) + 'x_' + str(j+1)
    print '方程的通解是：\n', 
    for i in range(n-1):
        print 'x_' + str(i+1), '=', result[i]            


# 得到简化的阶梯矩阵和主元列
def solve(a,b):
    sigma = arguemented_mat(a,b)
    print '增广矩阵为：'
    print sigma
    sigma, main_factor = trape_mat(sigma)
    sigma = back_solve(sigma, main_factor)
    print '方程的简化阶梯矩阵:'
    print sigma
    print '方程的主元列为：'
    print main_factor
    print_result(sigma, main_factor)
    return sigma, main_factor


if __name__ == '__main__':
    a = np.array([[1,6,2,-5,-2], [0,0,2,-8,-1], [0,0,0,0,1]])
    b = np.array([-4, 3, 7])
    sigma,mf  = solve(a,b)
    print '*' * 20
    a = np.array([[1,0,-5], [0,1,1], [0,0,0]])
    b = np.array([1,4,0])
    sigma,mf  = solve(a,b)
    print '*' * 20

增广矩阵为：
[[ 1  6  2 -5 -2 -4]
 [ 0  0  2 -8 -1  3]
 [ 0  0  0  0  1  7]]
方程的简化阶梯矩阵:
[[ 1  6  0  3  0  0]
 [ 0  0  1 -4  0  5]
 [ 0  0  0  0  1  7]]
方程的主元列为：
[0, 2, 4]
方程的通解是：
x_1 = 0-6x_2-3x_4
x_2 = x_2(free var)
x_3 = 5+4x_4
x_4 = x_4(free var)
x_5 = 7
********************
增广矩阵为：
[[ 1  0 -5  1]
 [ 0  1  1  4]
 [ 0  0  0  0]]
方程的简化阶梯矩阵:
[[ 1  0 -5  1]
 [ 0  1  1  4]
 [ 0  0  0  0]]
方程的主元列为：
[0, 1]
方程的通解是：
x_1 = 1+5x_3
x_2 = 4-1x_3
x_3 = x_3(free var)
********************