使用高斯消元法求解线性方程组

作为矩阵运算的一个应用，高斯消元法在线性方程组求解中发挥着重要的作用。高斯消元法的本质，是通过一系列线性变换，将线性方程组的系数矩阵变为上三角矩阵或者下三角矩阵，进而通过回代的方法得到所有变量的值。本文作为一篇抛砖引玉的介绍，将简要说一下如何在一个具体的线性方程组上应用高斯消元法，并给出一个python程序用于展示编程求解线性方程组的方法。

一、矩阵、矩阵的线性变换

矩阵的本质是一个数据表。例如，下面的这个上三角矩阵（主对角线以下都是零的矩阵）就是一个3×3的数据表。 $$\left[\begin{array}{ccc}1& 2& 3\\ 0& 5& 6\\ 0& 0& 9\end{array}\right]$$ 矩阵的线性变化包括初等行变换和初等列变换。初等行变换包括三种变换：（1）对换两行（例如，对换矩阵中$i$和$j$两行，记作$r_i\leftrightarrow r_j$）；（2）以一个非零实数乘某一行的所有元素（例如，第$i$行乘$k$，记作$r_i\times k$）；（3）把某一行的所有元素的$k$倍加到另一行上（例如第$j$行的$k$倍加到第$i$行，记作$r_i+k{r}_j$）。列变换也包括三种变化，将上述行变换的“行”换成“列”，就得到了初等列变换。

二、高斯消元法

高斯消元法(Gaussian elimination)是求解线性方阵组的一种算法。它通过逐步消除未知数来将原始线性系统转化为另一个更简单的等价的系统。它的实质是通过初等行变化，将线性方程组的增广矩阵转化为行阶梯矩阵。 我们举一个例子： $$\begin{array}{}\text{(1)}& 2x+y-z& =8\text{(2)}& -3x-y+2z& =-11\text{(3)}& -2x+y+2z& =-3\end{array}$$ 我们构造增广矩阵，也就是上述线性方程组的系数矩阵$A$加上常数向量$b$。 $$\left[\begin{array}{cccc}2& 1& -1& 8\\ -3& -1& 2& -11\\ -2& 1& 2& -3\end{array}\right]$$ 经过初等变换$r_2-(-3/2)r_1,r_3-(-1)r_2$得到 $$\left[\begin{array}{cccc}2& 1& -1& 8\\ 0& 1/2& 1/2& 1\\ 0& 2& 1& 5\end{array}\right]$$ 再经初等变换$r_3-4r_2$得到 $$\left[\begin{array}{cccc}2& 1& -1& 8\\ 0& 1/2& 1/2& 1\\ 0& 0& -1& 1\end{array}\right]$$ 于是我们得到了一个简化的三角方程组。在这个三角方程组中，首先就可以得到变量z的值，进而得到y的值，最后得到x的值。 $$\begin{gathered} \left(\begin{array}{}\text{(4)}& 2x+y-z& =8\text{(5)}& 1/2y+1/2z& =1\text{(6)}& -z& =1\end{array}\right) \\ \Rightarrow z=-1,y=3,x=2 \end{gathered}$$

三、程序实现

这里我们使用Python实现一个线性方程组的求解程序。这个程序的思想来自文章《高斯消元法详解》，并结合Python的一些特性进行了修改。

#!/usr/bin/python3
#coding=utf-8

import numpy as np

## ======== 以下为高斯消元法的算法部分 ========
# 参考代码：https://blog.csdn.net/lzyws739307453/article/details/89816311

# 打印矩阵的内容
def printM(a):
    m = len(a)
    n = len(a[0])
    for i in range(0,m):
        for j in range(0,n):
            print("%10f,\t"%a[i][j],end="")
        print()

# 选择列主元并消元
def SelectColE(a): #a:2-dimension matrix; n:int 
    n = len(a) # 参数n应该代表的是矩阵的行数，也就是线性方程组的方程个数
    for i in range(0,n): # 循环处理每一行的信息
        r = i
        for j in range(i,n):
            if(abs(a[j][i])>abs(a[r][i])): r=j
        if(r!=i): # 当r!=i时，交换两行的内容。这里的操作基于以下逻辑：
            # 由于单纯的高斯消元法会将矩阵化为梯形矩阵，靠下方的行中系数为0的项更多
            # 而非0项也会有很多分母，在除法运算中会导致精度降低
            # 因此将系数大的项尽可能的往前移动，可以提高精度
            tmp  = a[r]
            a[r] = a[i]
            a[i] = tmp
        for j in range(i+1,n): # 消元
            temp = a[j][i]/a[i][i]
            for k in range(i,n+1):
                a[j][k] -= a[i][k]*temp

# 高斯消元法
def Gauss(a): 
    n = len(a) # 参数n应该代表的是矩阵的行数，也就是线性方程组的方程个数
    SelectColE(a)
    for i in range(n-1,-1,-1): # 回代求解
        # i 从 n-1递减循环到0，可以取到0
        # 从第n-1行循环到第0行。由于前面已经对矩阵进行了变换，得到了一个倒三角矩阵
        # 所以n-1行可以直接求得第n-1个未知数的值（未知数的下标从0开始）
        # 而前面的行则需要经过一些迭代
        for j in range(i+1,n):  # 对行矩阵进行回代
            # 在此行的最后一个元素处存储回代后的值
            # a[i][n]是行矩阵的最后一个元素，是存储回代后的数值的地方
            # a[i][i]是本行对应的未知数的系数。一共i行，那么就是i个未知数
            # 例如一个3*4的矩阵，有3个变量分别是x0,x1,x2，那么当i=1时，
            # 对应位置就是变量x1的系数
            # 同样的，a[i][j]是本行第j个未知数的系数，
            # 而第j个未知数的值刚刚我们已经求过了，是a[j][n]，
            # 所以a[i][n]减去他俩相乘的值就行
            a[i][n] -= a[i][j]*a[j][n]
        # 回代结束后，对系数进行化简，使得第i个未知数的系数变为1，即可
        # 于是第i个未知数的值就等于a[i][n]/a[i][i]。我们把它存储在a[i][n]的位置上
        a[i][n] /= a[i][i]

## ======== 以上为高斯消元法的算法部分 ========

# 传入一个线性方程组的增广矩阵，这个函数将会调用高斯消元法的求解程序，并打印出求解结果。
def run(A):
    # 下面的内容是高斯消元法求解线性方程组的测试
    # 最终的方程组的解为a[:,-1]
    a = A.copy()
    print('a=\n',np.array(a))
    Gauss(a.copy())
    print("求解结果：")
    printM(a)
    print("数据结构：")
    printM(a)
    for i in range(0,len(a)):
        print("X%d = %9f"%(i,a[i][3]))
    print(np.array(a)[:,-1])
    print('\n\n')

# 主函数。此处展示了4个线性方程组的求解过程
if(__name__=='__main__'):
    print("Gauss_Elimination_Method.py")
    A1 =[[ 1,-1,-1, 2 ],
        [ 2,-1,-3, 1 ],
        [ 3, 2,-5, 0 ]]
    run(A1)
    
    A2= [[ 1, 2, 1, 8],
         [ 4, 0, 0,16],
         [ 0, 4, 0,12]]
    run(A2)

    A3= [[ 1, 2, 0, 8],
         [ 4, 0, 1,16],
         [ 0, 4, 0,12]]
    run(A3)

    A4= [[ 2, 0, 0, 8],
         [ 0, 1, 0,16],
         [ 4, 0, 1,12]]
    run(A4)