代数方程组求解｜高斯消去法
本文描述高斯消去法的基本处理思路。
有限体积法离散得到的线性方程组可以表示为下面的矩阵方程形式：
式中 为单元系数 构成的矩阵； 为未知量 构成的向量； 为源项 构成的向量。
式(1)矩阵形式为：
一般来说，矩阵中的每一行都表示计算域中某一单元上定义的局部离散方程，而其中的非零系数则与该单元的相邻单元相关。系数 可以衡量当前控制体形心处的 值预期相邻单元上的的关联程度。由于一个单元只与有限的几个单元相邻，其个数取决于离散域中单元的连接性，这便使得系数矩阵中很多元素都是零，因此该系数矩阵A总是稀疏的（即非零元素仅占矩阵中很小的一部分)。此外，如果使用结构网格系统，矩阵A将具有带状结构，即所有非零元素均出现在少量的几条对角线上。
1 引例
引入一个简单的例子。
考虑一个具有两个未知量 和 的二元线性方程组：
该方程组可以通过从其中一个方程中消去一个变量来实现求解，比如用式（4）去式（3）项乘以 可以消去 ，得到新的方程:
该方程只包含一个未知量 。通过式（5）得 ：
将得到的 回代至式（3）以求得 ：
上面的过程由两个步骤组成：
对方程进行操作消去一个未知数。这一步的最终结果是得到一个只含有一个未知量的方程（即一个方程一个未知量)。 直接对新方程求解，并将计算结果代回至其中一个方程中进而求得剩余的未知量。
同样的过程可以推广到由式（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》及《计算流体力学中的有限体积法》。
”
（完）
