将控制方程沿空间域和时间域分别离散后，可得到如下形式的代数方程

式中，a_P、a_N为系数，R_P为源项。

把求解域中所有单元的如式（11-10）的控制方程组装在一起，就构成了下面的代数方程组

AΦ=R（11-11）

式中，A是一个系数矩阵，系数a_P位于矩阵A的对角线上，系数a_N则位于对角线外；Φ是由所有单元的变量ϕ构成的张量；R为源项张量。

矩阵的稀疏程度取决于控制体积的排列形式，对角线上下的系数分别对应于单元面上的变量值。组装后的代数方程组（11-11）的求解方法主要有直接求解和迭代求解两种[267-269]。

1.代数方程组的直接解法

直接解法通过一系列的数学运算求得代数方程组的精确解。常用的直接解法有Gauss消去法、主元素消去法、Jordan消去法、Doolittle分解法等。直接解法一般适合于阶数不太高的线性方程组，其计算量较少，精度较高，但是程序较复杂。采用直接解法获得定解所需要的计算次数与代数方程组个数（或未知量个数）的立方成正比。这样，采用直接法获得最终解所需要的运算步骤随方程个数的增加而急剧增大，致使求解变得越来越昂贵和困难[267，268，270]。在实际应用中，方程组（11-11）是一个数量庞大的代数方程组，并且具有很强的非线性，因此在实际的模拟分析中，采用直接法求解几乎是不可能的。

2.代数方程组的迭代解法

与直接解法不同，迭代解法是用某种极限过程逐步逼近精确解的方法。在迭代法中，首先要假定一组初始值，然后根据方程逐步修正，直到解的误差减小到容差以内为止。迭代法通常能够识别和利用代数方程组系数矩阵的稀疏性进行计算，需要计算机的存储单元较少，程序相对简单。因此，在大型计算中，一般采用迭代解法。

以方程组（11-11）为例，假设迭代n次时方程组的近似解为Φⁿ，Φⁿ并不能完全满足方程组（11-11），而是有一个余量rⁿ，即

AΦⁿ=R_-rⁿ（11-12）

定义第n次迭代时解的迭代误差εⁿ为精确解Φ与近似解Φⁿ的差

εⁿ=Φ-Φⁿ（11-13）

用方程组（11-11）减去方程组（11-12），可知误差εⁿ满足下面的方程(www.daowen.com)

Aεⁿ=rⁿ（11-14）

方程组迭代求解的目的是使余量rⁿ逐渐减小，同时解的误差εⁿ逐渐减小到容差范围以内（或者减小到零）。假定迭代公式如下

MΦⁿ⁺¹=NΦⁿ+B（11_-15）式中，M、N和B为迭代公式的系数。当迭代第n步方程组的解收敛时，有Φⁿ⁺¹=Φⁿ=Φ，且解满足方程式（11-11），则有

（M-N）Φ=B （11-16）

并且

A=M-N

B=R （11-17）

或者，更一般地，各系数间满足下面的关系

PA=M-N

B=PR （11-18）

式中，P是一个非奇异矩阵，用来调节各系数间的关系。

一个有效的迭代方程式（11-15）需要满足计算容易、收敛快等要求。要使迭代计算容易，就要求NΦⁿ计算容易且系统易于求解。系数矩阵A和N都是稀疏矩阵，NΦⁿ的求解相对简单。要使迭代方程易于求解，就要求迭代矩阵M转置容易，通常迭代矩阵M为对角矩阵、三对角矩阵、三角形矩阵或五对角矩阵等。

最基本的迭代解法有Jacobi迭代法和Gauss-Seidel迭代法。这两种迭代方法均可非常容易地在计算机上实现，但是当方程组规模较大时，这两种方法收敛速度较慢^[271]。在本书后面的模拟分析中，主要采用ICCG迭代法^[272]来求解对称矩阵。对于非对称矩阵则用Bi-CG-STAB迭代方法来处理，其实现过程在Van Der Vorst的论文^[273]中给出了详细介绍。