2.4　椭圆方程和稳定状态问题

考虑椭圆式（2.3.18），利用中央差分方案进行求导，结果得到一个二阶的精度方案为

将式（2.3.18）降低成一个代数方程组为

通过替换I=j+（i-1）J；I=1，（I，J），该方程组可以转换为一个单一矩阵方程为

其中A是N×N的系数矩阵，X是一个长度为N的向量，它包含未知变量φ从i=1，j=1到i=I，j=J所有N=IJ个格网点处的值，B是相同长度的向量，它包含式（2.4.2）右侧的应力项。矩阵A是大而稀疏的，主要是有边缘的三对角矩阵。

许多CFD问题最终都降低为矩阵方程的解，例如式（2.4.3）。因此一种有效的“求逆”方法显得非常重要，尤其是矩阵非常大的时候（对于一个包含100×100格网的问题来说，矩阵尺寸为104×104），然而这种矩阵通常是稀疏矩阵，它的求解与拉普拉斯方程求解一样，因此可以对其进行简化。

对代数式（2.4.3）求解的两种主要方法是直接法和迭代法。直接法例如Gauss消元法和Gauss-Jordan消元法是计算密集型的，矩阵大而密集时，舍入误差也大，它们很少被用于大的三维CFD问题中。在矩阵稀疏且为带状矩阵的情形中，广泛使用的方法是特别版的Gauss消元法。但是大多数情况下，所选择的方法通常是迭代法，这种方法从一个初始估计开始，使用有效的加速手段计算到真实解。

2.4.1　直接解法

对式（2.4.3）直接求解的一种可用方法是传统的Gauss消元法或Gauss-Jordan消元法，这种消元过程使用系统标准化和系统减法将矩阵的所有子元素和超对角元素消除。可以对子对角元素进行高斯前向消元然后使用回带法，Press等（1992）对这些内容进行了讨论。消元过程要求每一步除以一个元素，称为枢轴点。为了使该方法尽可能精确，必须选择最大可能元素作为枢轴点，这通常要求在每一步中对行和列重排。将行进行交换后使最大元素转换到期望位置并不会改变最终的解，而将列交换则会使解变得杂乱。行交换称为部分轴元法，它执行起来较为简单，而列交换称为全部轴元法，它需要额外的计算，因为必须对杂乱的解重新组合。通常使用的方法是将原始方程标准化并使其最大系数为1的情况下选择可能是最大的元素作为枢轴点，这种方法称为隐式旋转。然而，这些方法要求方程右侧的向量B为已知。对全部轴元法的Gauss-Jordan消元法程序可以参考Press等（1992）。

2.4.1.1　LU分解法

LU分解法为直接法中的首选，由于它所得三角集的解的简单性质：替换。不论矩阵A是什么结构，无论它是稀疏的（大量元素为0）或者密集的（极少数或者没有0元素）或者是带状的（只有对角线上或其附近的元素为非零值），LU分解都要求矩阵A被分解成下三角阵L和上三角阵U：A=L·U，那么A·X=（L·U）X=L（U·X）=B，因此可以首先用简单的前向替换对L·Y=B进行求解，使得从L的结构可以找到向量Y，然后用Y求解U·X=Y，使可以用简单的后向替换从U的结构中得到X，也就是求解U·X=L-1B，之后的任务就是有效地进行LU分解。Press等（1992）列出了这样的算法，令

要得到系数αij和βij，写出N+N2方程如下

通过采用特殊顺序重组很容易对这些方程进行求解，使克劳特算法（Press等，1992）为

对每一个j=1，2，…，N ：

计算下一个j

绕轴旋转就是将一个矩阵的行和列重新排列，以选出最大可能的矩阵元素对矩阵的其他元素进行标准化，这对克劳特算法很重要。然而，只需对行进行交换的部分轴元法就已经足够，因为在矩阵求逆的过程中，对矩阵的行进行交换不会对解产生影响，它执行起来比全部轴元法简单，全部轴元法会使解变得混乱［对于该问题的讨论可以参考Press等（1992）］。如果L和U的元素事先已知，可以用前向替换对向量Y进行求解如下

然后可以用后向替换求未知的向量X

Press等（1992）和Fletcher（1998a）都列出了进行LU分解的FORTRAN源代码（分别为LUDCMP和FACT），还列出了前向替换和后向替换（分别为LUBKSB和SOLVE）。Press等使用更为高效的克劳特算法进行LU分解。这些实现都是对密集矩阵而言的，它们对于稀疏矩阵来说在存储空间和计算时间上都存在浪费问题，因为零元素也会被存储和计算。有专门为稀疏矩阵而设计的方法，它使得算术运算量和存储需求最小化。Press等（1992）为三对角矩阵提出了一种特别简单的最优化方法（TRIDAG），该方法是一种将分解和替换相结合的方法（该方法大部分情况下都不必使用绕轴旋转）。

LU分解的一个优势是，它不要求方程右侧的向量B为已知，因此对矩阵进行分解后，替换代码可以利用它来有效地对不同B向量进行重复求解，这也意味着对于很大的方程组来说，由舍入误差积累引起的严重精度丧失可以被修正，迭代改进对机器精度进行存储（Press等，1992）。令Xc为修正解，Ww为误差解，δX为误差。因为A·Xc=B，乘以A和Xw后得不到B，而是得到B+δB。因此，只通过替换就可以对A·δX=δB=（A·Xw-B）进行求解，从而得到修正的δX，而不需要再次对矩阵A进行分解，Press等（1992）为此提出了一个算法。

2.4.1.2　 Thomas算法

对单一未知变量的系统使用三点有限差分方法或者线性内插的有限元素法，结果得到具有三对角结构的矩阵，用分步法（本章2.5节）对多维问题进行求解也会产生三对角系统。线松弛迭代法也要求三对角矩阵的解（本章2.4.2节），因此三对角矩阵的有效算法对CFD很重要，一个特别有效的（但由于其具有递归性，矢量化和并行化比较困难）算法是Thomas算法，该算法由两个步骤组成：①是前向搜索，用矩阵主对角线上的单位元素将矩阵转化为上三角形式（在这个过程中向量B被修改为B′）②后向搜索，用简单回带法得到未知向量X，优化的Gauss消元法必须使用递归关系进行求解。

令aixi-1+bixi+cixi+1=di为三对角矩阵方程。进行扩展后得

矩阵A是非常稀疏的，只有主对角、次对角和超对角线上有值填充，不需要对整个N×N矩阵进行存储，只需要存储ai、bi和ci三个长度向量。一般的LU分解需要进行计算的次数是O（N）次，而不是O（N3/3）次。

Thomas算法使用的回归关系式为

其中系数αi和βi在前向搜索的过程中确定。将式（2.4.9）中的i改为i+1并将其替换为第i个方程，aixi-1+bixi+cixi+1=di，进行重新排列后得到

从而得到

这些是递归关系，αi和βi的值（可以由第一个方程b1x1+c1x2=d1进行确定）为

由此完成了前向搜索，递归关系式（2.4.9）中的所有系数都是已知的。从元素N进行的后向搜索提供了所要求的解，将式（2.4.9）带入到最后一个方程bNxN-1+bNxN=dN，得到

从式（2.4.14）中得到xN后，当i=N-1，N-2，…，1时，可以用式（2.4.9）确定xi，因为式（2.4.11）到式（2.4.12）中的系数都是已知的。

和所有的Gauss消元法一样，矩阵的病态条件会引起严重的舍入误差，为了避免这种情况发生，必须使|bi|＞|ai|+|ci|。对Thomas算法进行编码非常简单，Fletcher（1988a）列出了一个事例代码。

对包含单一未知变量的系统使用高阶精度方案时，会生成一个矩阵方程，该矩阵的带宽比三对角矩阵的带宽大。然而，Thomas算法可以被泛化从而处理这些仍然具有稀疏性的矩阵。Fletcher（1988a）对五对角矩阵的泛化Thomas算法进行T描述和讨论，如果未知变量多于一个，方程不止一个，那么求解过程通常会产生一个分块三对角矩阵，是M×M的子矩阵形成了三对角矩阵的对角线、次对角线和超对角线元素，而不是由标量形成。也可以将Thomas算法泛化到这种情形（Fletcher，1988a），它要求进行O（NM3）次运算。

在简单均匀矩形格网上对椭圆方程进行求解，矩阵元素是常数，在这种情形中，直接解法采用一个称为循环约化的过程，然后用傅里叶级数表示，这种方法是很高效的，因为它使用了现代的快速傅里叶转换（FFT）方法，不过这不是本书的讨论范围，这些方法的讨论读者可以参考Fletcher（1988a）和Press等（1992）。

2.4.2　迭代解法和松弛法

对于CFD和海洋模拟问题中的巨大稀疏矩阵特征来说，最常用的方法是迭代方法或松弛法。它将式（2.4.3）的矩阵A分裂成两部分，其中一部分很容易求逆（E），另一部分为残余部分（F）：A=E+F。那么方程变为

松弛法包括对如下方式进行连续迭代

相当于

迭代从初始估计X1。Rn是第n个迭代步骤的残差向量，当解收敛时它的范数趋于0。M=E-1F称为迭代矩阵，收敛速度取决于该矩阵的性质，尤其是它的谱半径σ，它是矩阵的最大特征值λi的大小。特征值是代数方程的根，通过对行列式进行归零处理得到

其中I是单位矩阵。然后对实数对称矩阵，所有的特征值都是实数。

谱半径的大小M决定收敛特性。如果σ＜1，则可以保证收敛性，这可以用误差en=Xt-Xn进行证明，该误差由方程en+1=Men=Mn·e0决定，因此在初始估计值e0处误差接近0，如果Mn接近0，则解是收敛的，如果谱半径小于1，则这一点可以保证，该条件通常转变为对矩阵A为对角占优矩阵的要求。谱半径的值越小，收敛速度越快。

2.4.2.1　Jacobi和Gauss-Seidel法

为了达到说明的目的，对拉普拉斯方程进行考虑，简单起见，假设x方向和y方向上的格网大小相同，那么式（2.4.2）变为

经典的点Jacobi（或Richardson）方法由迭代方案组成如下：

如果对应的矩阵A可分解为一个严格的下三角矩阵L（只有对角线下方的元素非零）、一个对角矩阵D（只有对角线上的元素非零）和一个严格的上三角矩阵U（只有对角线上方的元素非零），那么Jacobi方法等价于

或者等价于

对于Dirichlet边界条件，迭代矩阵M=D-1（L+U）的谱半径通常情况下格网大小不相等，则有

其中Δx和Δy是格网大小。通常情况下该值虽小于1但很接近于1，因此，实际上Jacobi方法的收敛速度较低。当I，J＞＞1时，假设Δx=Δy，σJ～1-（π2/2N），其中N=IJ为格网点的数量。收敛速度定义为将误差降低一个数量级所需的迭代次数，为-（logσJ）-1，因此它约为0.5N。因为每次迭代需要5N次运算，要使误差减少一个数量级需要2.5N2次运算。

加权Jacobi方法是式（2.4.20）的简单变体，即

其中可以选择ω使收敛在一定程度上得到加速。

点Gauss-Seidel（or Liebmann）法是Jacobi方法的一个简单变形如下

等价于

或

迭代矩阵M=（L+D）-1U的谱半径等于Jacobi谱半径的平方，因此，其收敛速度要快两倍，但是仍然较慢而不实用。Jacobi方法和Gauss-Seidel方法，收敛性都可以保证。

松弛法可以是上述的点松弛法，其中在每次迭代循环中对一个格网点进行更新。利用Gauss-Seidel松弛法作为例子，在式（2.4.25）的左侧唯一未知的点是（i，j），因此这是一种点迭代方法。假如认为格网点分布于一个棋盘上，很容易看出红色（偶数）点处的解只取决于黑色（奇数）点，反之亦然。因此，红色点集和黑色点集的交替更新是可能的，在每次松弛循环中的每次半搜索中对偶数点或奇数点进行更新，这一红黑点松弛法（或奇偶排序法）是一种最为常用的方法。

松弛法也可是线性松弛，其中每次对整行或整列的点进行松弛。对式（2.4.25）稍作修改得到

对每一个j处的所有i都可以获得一个解，也就是对j列的每一行所得三对角矩阵逐行求解，同样地，也可以进行逐列求解。这就是线Gauss-Seidel松弛法，它比点Gauss-Seidel法速度快两倍，收敛速度得到很大改善，而这是通过在每一步中引入点的解的信息，而不是相邻点的信息。

2.4.2.2　逐次超松弛（SOR）法

如果认识到这些松弛方法的本质是对之前的解添加一个校正这一事实的话，就很容易理解逐次超松弛（SOR）法的工作原理。SOR一直是进行椭圆求解的主要方法，直到最近才被现代多格网和共轭梯度法所取代。当使用点Gauss-Seidel方法时SOR的优势便显露出来，每步迭代中不是对校正项，而是对下方的式子进行

这等价于

或者

这是Gauss-Seidel值和之前的值的加权平均。α是松弛参数，收敛性要求0＜α＜2。对于超松弛法，要使它比Gauss-Seidel方法收敛快，则需满足α＞1。松弛参数的最优选择是此时谱半径为（Press等，1992）。需要注意，这仅是渐进有效的，只有在O（J）次迭代后得到。在初始阶段，误差逐渐增加。对于Jacobi和Gauss-Seidel方法采用的Dirichlet例，当Δx=Δy，I=J=N1/2时，将误差降低一个数量级所需要的迭代次数约为0.33N1/2。显然，该参数的选择很关键，因为在该最优值周围的一个狭窄窗口内收敛速度达到最好，离开该窗口后立即变差。不过，对一般问题而言，很难事先决定最优值，因此进行一些实验是很必要的。点松弛和线松弛，Jacobi和Gauss-Seidel松弛都可以是超松弛的，但对点Jacobi方法来说，松弛参数的最优值是1，也就是常规Jacobi方法已经是最优的。

每次迭代中α可以选择一个不同的值，因为上述最优值在开始时不是最优的。当α保持不变时，SOR方法称为静态超松弛法，而当α在每次迭代中不断改变时，称为非静态超松弛法。在迭代过程中使用一个特定的α序列使收敛速度得到极大增加，这种方法称为Chebyshev加速方法，Press等（1992）对该方法进行了描述，他在每次半搜索中改变α以保证误差指标单调减小。α值序列从1开始，以最优值结束：SOR的轻微改变形式称为对称SOR，式（2.4.17）中的矩阵E为（αL+D）D-1（U+D）/（2α-α2），它允许SOR可以被Chebyshev加速法进行加速（Fletcher，1988a）。SOR也可以是点SOR或线SOR，在最优松弛参数相同的情况下后者比前者快倍（Hirsch，1988）。

一个重要的问题是缓和，在该过程中可以在高性能计算机上对松弛方案进行向量化和并行化处理，比如常规的线Jacobi松弛法的向量化比线Ganss-Seidel SOR方法的向量化简单，因此，尽管它效率稍差，但也可能综合效果会更好。就像点松弛方法中的红黑方案那样，可以在连续搜索过程中在交替的线上执行Jacobi线松弛法，这称为斑马线松弛法，可以对它进行逐行或逐列执行，斑马线松弛法非常适合于现代计算机。

2.4.3　预条件共轭E梯度法

当使用SOR时，渴望使用一些方法使收敛得到加速。如前所述，收敛的速度取决于迭代矩阵的特征值λ，每个特征值表示一个因子，它使一次迭代中误差的特定本征模减小。显然，模的最大特征值决定收敛的速度。比值|λmax|/|λmin|是矩阵的条件数，表示特征值谱的带宽，如果它比1大得多，则特征值接近λmax占主导并降低收敛速度，对收敛进行改进的一种方法是用一个预处理矩阵乘以预条件矩阵，从而预条件矩阵的特征值谱的条件数较小，而在λ为λmax时最小。

另一种收敛加速方法是共轭梯度法，它要对如下的函数进行最小化

相当于令其梯度等于0

它的求解与式（2.4.3）一样。因此，方法简化为对n维空间中的函数进行最小化，这种方法称为共轭梯度法，它只适用于对称的正定矩阵（对所有的z来说，zTAz＞0），这是对椭圆方程求解的一般情形。这种方法由以下步骤组成（Fletcher，1988a）

其中Pn组成一个搜索方向（开始时等于初始半径R1），式（2.4.34）中的σn+1要保证这些搜索方向是相互正交的，而λn要保证式（2.4.32）在Pn方向上最小。在每一步中，基于A的特征函数，该方法通常使半径R的级数展开式的特征向量贡献消除。在没有舍入误差的情况下，该方法得到一个计算次数不超过N的鲁棒算法，如果A的一些特征值聚集在一起，它的计算次数更少，不过，舍入误差通常会降低这种方法的有效性。

如果对方程进行修改，使矩阵A可以转换为一个近似单位矩阵，就可以避免上述问题，这称为预处理，预条件共轭梯度法是从偏微分方程（PDEs）得到的有限差分方程进行求解的极有效方法（Fletcher，1998a；Press等，1992）。预处理减少了特征值的传播，因此收敛速度更快，这种方法等价于求解下列方程，而不是解式（2.4.3）

其中

选择矩阵W使WTW=E，选择不同的E会使W产生不同的方案，那么预条件共轭梯度方案归结为

其中S和代表伪残差。与之前一样，用σn+1的值表示收敛。

2.4.4　多格网法

使用现代多格网法对解加速的速度可能比使用SOR方法还快。可以认为，松弛法本质上是迭代地使不正确的初始估计所引起的误差场得到衰减，那么收敛的速度就取决于松弛方法使误差衰减的确切程度，松弛法在衰减误差场的高波数分量方面非常有效，通常只需要几次迭代就能使波数分量衰减到接近格网分辨率波数（2Δx）-1。不过，低波数的衰减非常缓慢，因此在几次迭代以后，总误差的进一步减小变得缓慢。多格网方法是通过对一个事实的巧妙探索而实现对松弛法的收敛进行加速的，这个事实就是误差的高频分量在松弛过程中被迅速消除，某个分辨率格网下的低频误差在另一较低分辨率格网下会表现为高频误差。取代在一个单一格网上进行求解，它们在分辨率为各格网分辨率乘积的多格网上成功地进行了求解。

可以使用一个简单的一维边界值问题对多格网方法进行说明（Briggs，1987）为

利用中央二阶精度方案，对应的有限差分方程变为

式（2.4.39）显然是内部格网点j=2，…，n-1处的解的一个（n-2）×（n-2）三对角系统，可以用算法对其进行有效求解，例如Thomas算法。更一般地说，多维椭圆问题的一个有限差分近似中得到的线性方程系统的形式如式（2.4.3），如果yj表示一个近似解，那么误差ej=xj-yj由残差方程A·E=R=B-A·Y控制，可以使B=0，然后得到精确解为X=0，误差就是E=-Y。

考虑y的初始估计由m=1到n的不同波数构成

用加权Jacobi方法进行100次迭代的应用表明，波数越高，衰减速度越快。波数分量n衰减得特别快，而波数1则衰减缓慢。而且，尽管前几次迭代中高波数分量的抑制使总误差大大减小，但误差的进一步降低却很缓慢，高频分量被迅速消除，而低频分量依然保留下来，因此，在前几次迭代后，收敛速度急剧降低。在整个迭代方案中，误差的低频分量消除缓慢使得收敛缓慢。松弛方法是有效的平滑器，然而，尽管它们能对初始误差分量进行有效抑制，却使平滑分量本质上保持不变。收敛的速度不会受到加权因子的太大影响，实际上反而会受格网空间h的任何减小的影响。因此，尝试重新定义格网以获得更精确的解是不可能的，反而只会起到反效果。

多格网方法的指导原则是，在松弛过程中，振荡分量（m＞＞n/2）大幅度衰减，而平滑分量（m＜＜n/2）则不会（Brandt，1977；Briggs，1987；Hackbusch，1985；McCormick，1987）。假设误差的振荡分量在某一特定格网中已经被利用几次松弛搜索进行消除。因为平滑分量相对保持不变，进一步迭代不会带来什么优势，但是对密格网来说已经平滑的分量对粗格网来说却是振荡的，在粗格网上进行松弛会使这些误差得到一定程度的消除，至少会变得更快。因此，用粗格网可能会更好，在该格网中的点上进行松弛，同时也进行松弛搜索。

这个过程可以从获得期望解的格网开始，在进行几次迭代后，当收敛变慢时，残差被转换到一个粗格网中。由于密格网上误差的低频分量在粗格网上看起来似乎是高频分量，在对这些误差分量进行衰减时，在粗格网上进行几次迭代比在密格网上进行更多次的连续迭代要有效得多，然后反过来将这一改进的残差估计（对残差方程进行松弛处理）转换到密格网中，从而获得更好的误差估计并使最终得到的解更精确，这一在密格网和粗格网上进行连续求解的过程可以保持收敛性，这种方法称为粗格网校正，该过程也可以持续到格网达到最粗糙的情况（2×2），或者也可在某些中间点处终止。

另一种方法是从粗格网开始进行松弛搜索，并将其解作为更密格网上的一个初始估计，接着又再作为下一个更密格网的初始估计，直到得到满意的解为止，重复整个循环直到收敛，这被称为嵌套迭代。两种方法都要求在两个格网间对信息进行正确且高效转换的方法，即从粗糙到密集和从密集到粗糙，前者称为延伸，后者称为限制。通常对延伸来说最合适的是简单的线性内插法，而对于限制来说，则较合适的是简单注入法，这种方法中粗格网点假设对应的密格网点会更好（尽管周围的密集格网点的加权平均也会很好）。一个简单的双重粗格网校正方案如下（一直重复直到收敛）：

（1）用初始估计Yh在密格网Ωh上对AhXh=Bh进行p次松弛。

（2）在粗格网Ω2h上计算残差

（3）用初始估计E2h=0在粗格网Ω2h上对A2hE2h=R2h进行p次松弛。

（4）修正密格网分辨率

（5）用新的初始估计Yh在密格网Ωh上对AhXh=B h进行p次松弛。

其中上标h表示原始格网，2h表示格网尺寸比原始格网大两倍的粗格网。这个双重格网方法是非常高效的，渐进收敛速度与格网大小有关。工作量与格网点总数N=nd大小相当，n为d问题维之一的格网点数，对最优的SOR来说为N3N2，对共轭梯度方法来说为N5/4

2.4.4.1　V循环

为了对粗格网（或分辨率更小的格网）平滑特性的优势进行探索，采用的格网可能大于2个：

用初始估计Yh在密格网Ωh上对AhXh=Bh进行p次松弛。

在粗格网Ω2h上计算残差

用初始估计E2h=0在粗格网Ω2h上对A2hE2h=R2h进行p次松弛。

在粗格网Ω4h上计算残差

用初始估计E4h=0在粗格网Ω4h上对A4h E4h=R4h进行p次松弛。

在粗格网Ω8h上计算残差

…

在粗格网ΩLh上对ALhELh=RLh进行松弛。

…

修正密格网分辨率

用初始估计Y4h在密格网Ω4h上对A4h X4h=B4h进行q次松弛。

修正密格网分辨率

用新的初始估计Y2h在密格网Ω2h上对A2h E2h=R2h进行q次松弛。

修正密格网分辨率

用新的初始估计Yh在密格网Ωh上对AhXh=B h进行q次松弛。

其中，上标4h表示格网大小为原始格网的四倍的粗格网，以此类推。这种方法称为V循环多格网方法，它非常高效，多格网方法有很多可能的变形，比如W循环，与效率最好的基于FFT的直接解法的效率O（ndlogn）相比，该方法的效率为O（nd）。然而，必须注意的是，SOR不能用在多格网方法中，因为对多格网方法很关键的高频衰减能力比常规松弛方法（例如Gauss-Seidel法）差很多，通常最适合多格网方法的是黑红点松弛。更多细节和软件方面的信息，读者可以参考Fletcher（1988a）和Press等（1992）的著作。注意，对复杂形状域（比如海盆）来说，基于FF T的方法和多格网方法的应用都很复杂。

需要注意，虽然松弛可应用到线性系统的原始方程或残差方程中，但如果系统是非线性的，则必须将它用到原始方程中，更多细节读者可以参考Brandt（1977）、Hackbush（1985）、Briggs（1987）、Fletcher（1988a）和Press（1992）等。

2.4.5　伪瞬态方法

稳态边界值问题迭代方法通常等价于推进一个对应的时间依赖问题到它的渐进状态。要证明这一点，可以将时间导数添加到拉普拉斯式（2.3.18）以得到一个双曲线初始值问题

该方程可以用一个任意初始状态求解到它的稳定状态，与原始边界值问题一样受相同边界条件约束，因此渐进稳定状态对应于原始边界值问题的解。使用显式FTCS（时间前向，空间中央）方案（并假设Δx=Δy），得到

式中：n表示当前时间，Δt是时间步长，要满足稳定性，（本章2.5节），当它等于0.25时，除了n表示时间级而不是迭代之外，所得到的方程与式（2.4.20）Jacobi迭代方案完全一样。

这种人工添加一个时间依赖项从而将边界值问题转换为初始值问题并求渐进稳态解的方法称为伪瞬态方法，它是一种克服求解复杂稳态问题的有效方法，复杂稳态问题中方程的性质有可能会在部分域中从双曲线到椭圆发生改变，该方法在20世纪60年代早期被成功应用于当时几乎不可处理的钝再入体周围高马赫数流的稳态问题，这种流在驻点附近加速回到超音速状态之前立即减速到次音速状态。因此，这里的域包含嵌入的椭圆区域，而激震前沿首先是未知的。不同域中跨越边界求解的复杂性必须作为解的一部分，这种复杂性可以避免，方法是通过将稳态问题转换为一个初始值问题并从估计的初始状态向前推进到渐进稳态。求解渐进状态的等价非稳定问题有一个很大优势，即方程是双曲线的，可以暂时使用多维时间依赖问题的强大分割技术，比如交替方向隐式技术（ADI）。