本文描述高斯消去法的基本处理思路。

有限体积法离散得到的线性方程组可以表示为下面的矩阵方程形式：

式中为单元系数构成的矩阵；为未知量构成的向量；为源项构成的向量。

式(1)矩阵形式为：

一般来说，矩阵中的每一行都表示计算域中某一单元上定义的局部离散方程，而其中的非零系数则与该单元的相邻单元相关。系数可以衡量当前控制体形心处的值预期相邻单元上的的关联程度。由于一个单元只与有限的几个单元相邻，其个数取决于离散域中单元的连接性，这便使得系数矩阵中很多元素都是零，因此该系数矩阵A总是稀疏的（即非零元素仅占矩阵中很小的一部分)。此外，如果使用结构网格系统，矩阵A将具有带状结构，即所有非零元素均出现在少量的几条对角线上。

1 引例

引入一个简单的例子。

考虑一个具有两个未知量和的二元线性方程组：

该方程组可以通过从其中一个方程中消去一个变量来实现求解，比如用式（4）去式（3）项乘以可以消去，得到新的方程:

该方程只包含一个未知量。通过式（5）得：

将得到的回代至式（3）以求得：

上面的过程由两个步骤组成：

1. 对方程进行操作消去一个未知数。这一步的最终结果是得到一个只含有一个未知量的方程（即一个方程一个未知量)。
2. 直接对新方程求解，并将计算结果代回至其中一个方程中进而求得剩余的未知量。

同样的过程可以推广到由式（1）式（2）所描述的N维方程组中。

2 消元

在接下来的推导中，的第1行表示的离散方程，第⒉行表示的方程，第行表示的方程。首先，从中第1行以下的所有方程中消去。为了能从第行中消去，需要用乘以第1行的系数，并用第行减去新得到的方程。执行完此步骤后，方程组变为：

然后在新的中把第⒉行以下的所有方程中的消去。为了将第行()中的消去，需要用乘以第2行的系数，并用第行减去新得到的方程。然后在新的系数矩阵中把第3行以下的所有方程中的消去，并重复执行此操作，知道将第行中的消去。至此，将得到一个等价方程组，如下所示，其矩阵被转化成一个上三角矩阵:

上述消元过程可以使用下面的算法进行描述：

for k = 1 to N-1
{
    for i = k+1 to N
    {
        ratio = a_ik/a_kk
        {
            for j = k+1 to N
                a_ij = a_ij - ratio * a_kj
        }
        b_i = b_i - ratio * b_k
    }
}

3 回代

从式（9）以看到，第个方程只有唯一未知量，因此利用该方程可以求得：

第个方程包含两个未知量和。既然已经得到，则可以利用得到的求，结果如下：

继续向回执行上述操作，直到到达第个方程，此时变量已经变为已知量，则的计算公式如下：

继续该操作，直到求得。

回代过程的算法总结为：

phi_N = b_N / a_NN
for i = N-1 to 1
{
    term = 0
    {
        for j = i+1 to N
        {
            term = term + a_ij * phi_j
        }
    }
    phi_i = (b_i - term)/a_ii
}

高斯消去法的计算量非常大，求解维线性方程组的运算量与成正比，其中的运算量用于回代过程。

注：对于使用合适离散算法得到的CFD离散方程系数矩阵，其对角线元素不为零，且矩阵满足对角占优，因此理论上讲都可以直接利用消去法求解。

”

高斯消去法编程示例（程序参考https://zhuanlan.zhihu.com/p/386954541）：

import numpy as np
 
def gauss(a, b):
    m, n = a.shape  # 获取矩阵的行数和列数
    c = np.zeros(n)  # 根据矩阵的行数构建一个一维0数组
    for i in range(n):
        # 限制条件
        if (a[i][i] == 0):  # 用高斯消去法解线性方程组时对角线元素不能为0
            print("no answer")
 
    # k表示第一层循环，(0，n-1)行
    # i表示第二层循环,(k+1,n)行,计算该行消元的系数
    # j表示列
    for k in range(n - 1):
        for i in range(k + 1, n):
            c[i] = a[i][k] / a[k][k]  # 计算出系数
 
            for j in range(k,m):  # 从K开始，减少不必要的计算
                a[i][j] = a[i][j] - c[i] * a[k][j]  # 对矩阵进行高斯消去
            b[i] = b[i] - c[i] * b[k]
 
        print('step',k,':n',a,'n')
        #print(b)
    x = np.zeros(n)
 
    x[n - 1] = b[n - 1] / a[n - 1][n - 1]  # 解出x[n-1],为回代作准备
 
    # 回代求出方程解
    for i in range(n-2, -1, -1):
        sum= 0.0
        for j in range(n-1, -1, -1):
            sum= sum + a[i][j] * x[j]
        x[i] = (b[i]-sum) / a[i][i]
    # 输出计算结果
    for i in range(n):
        print("x" + str(i + 1) + " = ","%.2f" % x[i])  # 输出结果
 
if __name__ == '__main__':
    a = np.array([
        [2.0, -1.0, 3.0, 2.0],
        [3.0, -3.0, 3.0, 2.0],
        [3.0, -1.0, -1.0, 2.0],
        [3.0, -1.0, 3.0, -1.0]
    ])
    b = np.array([6.0, 5.0, 3.0, 4.0])
    gauss(a, b)

程序输出结果：

step 0 :
 [[ 2.      -1.       3.           2. ]
 [ 0.      -1.5     -1.5     -1. ]
 [ 0.       0.5     -5.5     -1. ]
 [ 0.       0.5     -1.5     -4. ]] 

step 1 :
 [[ 2.             -1.               3.       2.        ]
 [ 0.             -1.5            -1.5            -1.        ]
 [ 0.              0.      -6.             -1.33333333]
 [ 0.              0.               -2.             -4.33333333]] 

step 2 :
 [[ 2.             -1.                  3.                  2.        ]
 [ 0.             -1.5              -1.5                -1.        ]
 [ 0.              0.                 -6.                 -1.33333333]
 [ 0.              0.                  0.                 -3.88888889]] 

x1 =  1.00
x2 =  1.00
x3 =  1.00
x4 =  1.00

可以看到每一步消去过程后的系数矩阵以及最终求解得到的结果。

高斯消去法处理逻辑比较简单，程序也比较容易编写。然而此方法随着方程组规模的增大，其计算量成指数方式增长。而且高斯消去法在求解过程中需要将整个系数矩阵读入内存，故对内存的需求也非常大，因此在实际的工程应用中不太可能使用此方法。

注：本文内容采集自《The Finite Volume Method in Computational Fluid Dynamics》及《计算流体力学中的有限体积法》。

”

（完）