FVM in CFD 学习笔记_第10章_补充专题_多重网格算法

参看资料：
https://ocw.mit.edu/courses/mathematics/18-086-mathematical-methods-for-engineers-ii-spring-2006/readings/am63.pdf
https://computing.llnl.gov/projects/hypre-scalable-linear-solvers-multigrid-methods/CiSE_2006_amg_220851.pdf
user.it.uu.se/~maya/Projects/3phase/AMG_parallel_Falgout.pdf

由于书中所讲的多重网格算法过于简略，无法参悟透彻，故重开一篇，再讲讲如何使用该算法。看了网上很多资料，大多语焉不详，要么是大部头的理论剖析，要么是寥寥数语不知如何应用，确实让人稀里糊涂无法弄得很透彻。本篇纯属作者个人胡蒙乱猜的，不当错误之处在所难免，还请及时指正，免得以讹传讹，贻误后人。闲言少叙，开吹：

在传统的迭代方法中，迭代是在固定网格上进行的，这样就会使得误差中的高频分量（波长较短的部分）衰减较快，而低频分量（波长较长）的部分衰减较慢，而多重网格的思想是让粗细不同网格层级之间来回计算折腾，从而使得各种高低频率长短波长的误差分量都得到快速地衰减，那么迭代收敛速度比传统的迭代方法要快得多了。

多重网格算法分为几何多重网格和代数多重网格，几何多重网格是真的划分了网格，真的构造了粗细网格层级；而代数多重网格则是直接跳脱出了网格的束缚，仅仅根据得出的系数矩阵$\mathbf{A}$，运用多重网格的思想，来进行虚拟的网格层级处理，从而提高收敛速度。

显然，几何多重网格更容易理解，也更简单一些，而在几何多重网格中，结构化的网格处理起来则更为方便一些，那么咱们就先从结构网格开始吧。

1 结构网格上的多重网格算法

由于要把粗网格和细网格的信息相互传递转化，所以实际上需要3个矩阵来做这些事情（用下标$h$来标识细网格，用下标$2h$来标识粗网格）：

• 把矢量从细网格到粗网格传递的限制矩阵${\mathbf{R}}_h^{2h}$
• 反过来，把矢量从粗网格向细网格延拓的插值矩阵${\mathbf{I}}_{2h}^h$
• 当然了，还要根据细网格的系数矩阵${\mathbf{A}}_h$，来构造粗网格上的系数矩阵${\mathbf{A}}_{2h}$，那么这个矩阵实际上就是${\mathbf{A}}_{2h}={\mathbf{R}}_h^{2h}{\mathbf{A}}_h{\mathbf{I}}_{2h}^h$

1.1 插值矩阵${\mathbf{I}}_{2h}^h$

先来看插值矩阵${\mathbf{I}}_{2h}^h$，对于1维问题，结构网格，计算域$0\le x\le 1$，边界值为0，用$h=1/8$来剖分细网格，用$2h=1/4$来剖分粗网格，边界点的值不参与编码和计算，则细网格上的变量$\mathbf{u}=[u_1,u_2,u_3,u_4,u_5,u_6,u_7{]}^T$共7个内部节点，而粗网格上的变量$\mathbf{v}=[v_1,v_2,v_3{]}^T$共3个内部节点。

如何根据粗网格（$2h$）上的变量$\mathbf{v}$来插值算得细网格（$h$）上的变量$\mathbf{u}$呢？不难发现$u_2,u_4,u_6$是和$v_1,v_2,v_3$重合的节点上的变量值，所以直接让它们相等就好了。而对于$u_1,u_3,u_5,u_7$这些点处的变量值，它们的位置是位于两个粗单元节点之间的（且是中点），所以用其两侧的粗单元节点值的算术平均来处理就好了（注意，边界节点的值是0），即$u_1=v_1/2$，$u_3=v_1/2+v_2/2$，$u_5=v_2/2+v_3/2$，$u_7=v_3/2$。写成矩阵形式，即
$$\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\\ u_6\\ u_7\end{array}\right]=\frac{1}{2}\left[\begin{array}{ccc}1& & \\ 2& & \\ 1& 1& \\ & 2& \\ & 1& 1\\ & & 2\\ & & 1\end{array}\right]\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]$$
从而获得了插值矩阵${\mathbf{I}}_{2h}^h$
$${\mathbf{I}}_{2h}^h=\frac{1}{2}\left[\begin{array}{ccc}1& & \\ 2& & \\ 1& 1& \\ & 2& \\ & 1& 1\\ & & 2\\ & & 1\end{array}\right]$$

1.2 限制矩阵${\mathbf{R}}_h^{2h}$

如果知道了细网格上的量，又如何将其限制到粗网格上呢？一种简单的想法是让相同节点的粗细网格值相互传递，即，$[v_1,v_2,v_3]=[u_2,u_4,u_6]$，这样做虽然看起来蛮好的，但问题在于，没有充分利用到细网格上面所有节点的信息，即，没有考虑到$[u_1,u_3,u_5,u_7]$节点值的影响。常用的处理方法是根据之前的插值矩阵${\mathbf{I}}_{2h}^h$来做反向处理，即让${\mathbf{R}}_h^{2h}=\frac{1}{2}({\mathbf{I}}_{2h}^h{)}^T$，称为全加权算子${\mathbf{R}}_h^{2h}$，如此一来
$$\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]=\frac{1}{4}\left[\begin{array}{ccccccc}1& 2& 1& & & & \\ & & 1& 2& 1& & \\ & & & & 1& 2& 1\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\\ u_6\\ u_7\end{array}\right]$$
写成分量形式
$$v_i=\frac{1}{4}[u_{2i-1}+2u_{2i}+u_{2i+1}],i=1,2,3$$
即每个粗网格节点的值，由与其重合的细网格节点值和其邻近的细网格节点值加权得到。全加权算子${\mathbf{R}}_h^{2h}$为
$${\mathbf{R}}_h^{2h}=\frac{1}{4}\left[\begin{array}{ccccccc}1& 2& 1& & & & \\ & & 1& 2& 1& & \\ & & & & 1& 2& 1\end{array}\right]$$

1.3 二维问题的插值矩阵和限制矩阵

还是结构网格，如果问题扩展为二维，即笛卡尔网格，那么这时候插值矩阵和限制矩阵是啥样的呢？

用$u_{i,j}$和$v_{i,j}$分别表示密网格和粗网格上的变量，先看从粗网格到细网格转化变量的插值矩阵，跟一维问题类似，如果细网格节点和粗网格节点重合，那么直接让两者存储的变量相同，即
$$u_{2i,2j}=v_{i,j}$$
如果细网格节点刚好落在粗网格两节点之间，就用粗网格这俩节点上变量平均来计算，即
$$\begin{gathered} u_{2i+1,2j}=\frac{1}{2}(v_{i,j}+v_{i+1,j})u_{2i-1,2j}=\frac{1}{2}(v_{i,j}+v_{i-1,j}) \\ {u}_{2i,2j+1}=\frac{1}{2}(v_{i,j}+v_{i,j+1})u_{2i,2j-1}=\frac{1}{2}(v_{i,j}+v_{i,j-1}) \end{gathered}$$
如果细网格节点是位于四个粗网格节点之间的，那么用这四个粗网格节点上变量的平均来计算，即
$$\begin{gathered} u_{2i+1,2j+1}=\frac{1}{4}(v_{i,j}+v_{i+1,j}+v_{i,j+1}+v_{i+1,j+1}) \\ {u}_{2i+1,2j-1}=\frac{1}{4}(v_{i,j}+v_{i+1,j}+v_{i,j-1}+v_{i+1,j-1}) \\ {u}_{2i-1,2j+1}=\frac{1}{4}(v_{i,j}+v_{i-1,j}+v_{i,j+1}+v_{i-1,j+1}) \\ {u}_{2i-1,2j-1}=\frac{1}{4}(v_{i,j}+v_{i-1,j}+v_{i,j-1}+v_{i-1,j-1}) \end{gathered}$$
如此，便可推得插值矩阵${\mathbf{I}}_{2h}^h$。

二维问题的限制矩阵${\mathbf{R}}_h^{2h}=\frac{1}{4}({\mathbf{I}}_{2h}^h{)}^T$，粗网格的每个节点值由细网格上的9个节点值插值得到，各节点系数如图所示。

即
$$v_{i,j}=\frac{1}{16}[4u_{2i,2j}+2(u_{2i+1,2j}+u_{2i-1,2j}+u_{2i,2j+1}+u_{2i,2j-1})+(u_{2i+1,2j+1}+u_{2i+1,2j-1}+u_{2i-1,2j+1}+u_{2i-1,2j-1})]$$
这样，也可推得限制矩阵${\mathbf{R}}_h^{2h}$。

对于三维问题，也可用类似方法推出插值矩阵和限制矩阵。

1.4 粗网格系数矩阵${\mathbf{A}}_{2}h$

粗网格上的系数矩阵${\mathbf{A}}_{2h}$是由细网格上的系数矩阵${\mathbf{A}}_h$和变量传递关系转化而来的，故有${\mathbf{A}}_{2h}={\mathbf{R}}_h^{2h}{\mathbf{A}}_h{\mathbf{I}}_{2h}^h$，完美。

比如，如果一维问题的控制方程和零边值条件为
$$\begin{gathered} -\frac{d^2\varphi }{d{x}^2}=f(x),0\le x\le 1 \\ \varphi (x=0)=\varphi (x=1)=0 \end{gathered}$$
跟前面讲的网格剖分方式一样，细网格等分8个单元，单元长度为$h=1/8$，用中心格式近似2阶项，有
$$-\frac{d^2\varphi }{d{x}^2}=-\frac{{\varphi }_{i+1}-2{\varphi }_i+{\varphi }_{i-1}}{h^2}$$
则细网格上系数矩阵${\mathbf{A}}_h$为
$${\mathbf{A}}_h=\frac{1}{h^2}\left[\begin{array}{ccccccc}2& -1& & & & & \\ -1& 2& -1& & & & \\ & -1& 2& -1& & & \\ & & -1& 2& -1& & \\ & & & -1& -2& -1& \\ & & & & -1& 2& -1\\ & & & & & -1& 2\end{array}\right]$$
粗网格等分4个单元，则粗网格上系数矩阵为
$$\begin{array}{rl}{\mathbf{A}}_{2h}& ={\mathbf{R}}_h^{2h}{\mathbf{A}}_h{\mathbf{I}}_{2h}^h\\ & =\frac{1}{4}\left[\begin{array}{ccccccc}1& 2& 1& & & & \\ & & 1& 2& 1& & \\ & & & & 1& 2& 1\end{array}\right]\frac{1}{h^2}\left[\begin{array}{ccccccc}2& -1& & & & & \\ -1& 2& -1& & & & \\ & -1& 2& -1& & & \\ & & -1& 2& -1& & \\ & & & -1& -2& -1& \\ & & & & -1& 2& -1\\ & & & & & -1& 2\end{array}\right]\frac{1}{2}\left[\begin{array}{ccc}1& & \\ 2& & \\ 1& 1& \\ & 2& \\ & 1& 1\\ & & 2\\ & & 1\end{array}\right]\\ & =\frac{1}{(2h{)}^2}\left[\begin{array}{ccc}2& -1& \\ -1& 2& -1\\ & -1& 2\end{array}\right]\end{array}$$

1.5 多重网格算法流程

如果仅有两层网格，细网格-粗网格-细网格这样做完一个循环，就是两层网格的V型循环，其算法流程如下：

1. 在细网格上迭代求解${\mathbf{A}}_h\mathbf{u}={\mathbf{b}}_h$，用Jacobi或Gauss-Seidel迭代，注意只做3次迭代步，得到尚未收敛的解${\mathbf{u}}_h$（虽然没收敛，但是高频分量衰减了很多）；
2. 计算细网格上面的残差${\mathbf{r}}_h={\mathbf{b}}_h-{\mathbf{A}}_h{\mathbf{u}}_h$，并将${\mathbf{r}}_h$传递到粗网格上面去，用${\mathbf{r}}_{2h}={\mathbf{R}}_h^{2h}{\mathbf{r}}_h$得到粗网格上面的残差${\mathbf{r}}_{2h}$；
3. 计算粗网格上系数矩阵${\mathbf{A}}_{2h}={\mathbf{R}}_h^{2h}{\mathbf{A}}_h{\mathbf{I}}_{2h}^h$（这个可以放到最前面算好存好拿来直接用，不用每次都算的），求解粗网格上方程${\mathbf{A}}_{2h}{\mathbf{E}}_{2h}={\mathbf{r}}_{2h}$，同样用迭代方法来求解，还是迭代3步就好，迭代初值选择为$\mathbf{E}=0$，迭代完成后算出了修正值${\mathbf{E}}_{2h}$；
4. 将粗网格上的修正值${\mathbf{E}}_{2h}$插值回细网格，用${\mathbf{E}}_h={\mathbf{I}}_{2h}^h{\mathbf{E}}_{2h}$，获得细网格上的修正值${\mathbf{E}}_h$；
5. 在细网格上继续求解${\mathbf{A}}_h\mathbf{u}={\mathbf{b}}_h$，还是用迭代方法迭代3步，注意这时候迭代初值选择为${\mathbf{u}}_h+{\mathbf{E}}_h$，即包含了粗网格的修正成分。

多重网格一般会构造好多层的粗细网格，而粗细网格遍历循环方式可以由V、W、F型，如图

1.6 例子

如果没有例子，那么还是搞不明白，这时候可以找个很简单的小例子来试试看，就用前面提到的一维零边值问题，控制方程和边界条件如下
$$\begin{gathered} -\frac{d^2\varphi }{d{x}^2}=f(x),0\le x\le 1 \\ \varphi (x=0)=\varphi (x=1)=0 \end{gathered}$$
用前面提到的网格剖分方式，细网格等分8单元，粗网格等分4单元，插值矩阵、限制矩阵、细网格系数矩阵、粗网格系数矩阵前面都推出来了。我们让方程中的$f(x)=-\sin (2\pi x)$，那么这个问题的精确解就很容易得到是$\varphi (x)=-\sin (2\pi x)/(2\pi {)}^2$，即细网格上的精确解
$$\mathbf{u}=[u_1,u_2,u_3,u_4,u_5,u_6,u_7{]}^T=[-0.0179,-0.0253,-0.0179,-0.0000,0.0179,0.0253,0.0179{]}^T$$
咱们用最简单的V型粗细两层多重网格算法来算算看，以Gauss-Seidel迭代为例（可以写个小程序来做哈，很简单的）。

1. 选择初始值$\mathbf{u}=0$，Gauss-Seidel迭代方法迭代3步，求解${\mathbf{A}}_h\mathbf{u}={\mathbf{b}}_h$，得到${\mathbf{u}}_h$为
   $${\mathbf{u}}_h=[-0.0122,-0.0172,-0.0122,-0.0000,0.0122,0.0172,0.0122{]}^T$$
2. 计算细网格上面的残差${\mathbf{r}}_h={\mathbf{b}}_h-{\mathbf{A}}_h{\mathbf{u}}_h$，得
   $${\mathbf{r}}_h=[-0.2500,-0.3536,-0.2500,-0.0000,0.2500,0.3536,0.2500{]}^T$$
   将${\mathbf{r}}_h$传递到粗网格上面去，用${\mathbf{r}}_{2h}={\mathbf{R}}_h^{2h}{\mathbf{r}}_h$得到粗网格上面的残差${\mathbf{r}}_{2h}$
   $${\mathbf{r}}_{2h}=[-0.3018,-0.0000,0.3018{]}^T$$
3. 求解粗网格上方程${\mathbf{A}}_{2h}{\mathbf{E}}_{2h}={\mathbf{r}}_{2h}$，Gauss-Seidel迭代方法迭代3步，注意迭代初值选择为$\mathbf{E}=0$，迭代完成后算出了修正值${\mathbf{E}}_{2h}$
   $${\mathbf{E}}_{2h}=[-0.0094,-0.0000,0.0094{]}^T$$
4. 将粗网格上的修正值${\mathbf{E}}_{2h}$插值回细网格，用${\mathbf{E}}_h={\mathbf{I}}_{2h}^h{\mathbf{E}}_{2h}$，获得细网格上的修正值${\mathbf{E}}_h$
   $${\mathbf{E}}_h=[-0.0047,-0.0094,-0.0047,-0.0000,0.0047,0.0094,0.0047{]}^T$$
5. 在细网格上继续求解${\mathbf{A}}_h\mathbf{u}={\mathbf{b}}_h$，注意这时候迭代初值选择为${\mathbf{u}}_h+{\mathbf{E}}_h$
   $${\mathbf{u}}_h={\mathbf{u}}_h+{\mathbf{E}}_h=[-0.0169,-0.0267,-0.0169,-0.0000,0.0169,0.0267,0.0169{]}^T$$
   显而易见，粗网格的修正起了作用，这时候的${\mathbf{u}}_h$更加接近精确解了，采用Gauss-Seidel迭代方法再迭代3步，得
   $${\mathbf{u}}_h=[-0.0189，-0.0257，-0.0189，-0.0000，0.0189，0.0257，0.0189{]}^T$$
   其残差${\mathbf{r}}_h={\mathbf{b}}_h-{\mathbf{A}}_h{\mathbf{u}}_h$的2范数为0.2165。还可以继续进行该细-粗-细的循环直到解收敛。
   注意，如果不用多重网格，只用Gauss-Seidel迭代方法只在细网格上进行迭代，进行了6步迭代后，所得到的的值${\mathbf{u}}_h$为
   $${\mathbf{u}}_h=[-0.0165，-0.0233，-0.0165，-0.0000，0.0165，0.0233，0.0165{]}^T$$
   其残差${\mathbf{r}}_h={\mathbf{b}}_h-{\mathbf{A}}_h{\mathbf{u}}_h$的2范数为0.2500，跟多重网格相比（0.2165）还是有差距的，表明多重网格方法的确是能起到加速收敛的作用。

2 代数多重网格（AMG）算法

把非结构网格上的多重网格算法，跟代数多重网格算法放到一起来讲，两者并没有区别。因为非结构网格的节点邻近拓扑关系是未知的、无规则的，而代数多重网格则更是直接刨掉了网格的束缚，然而两者的系数矩阵是相同的，即都是稀疏矩阵，哪里出现非零元素是无规则的，那么它俩放一块讲就是了。

注意，系数矩阵$\mathbf{A}$仍然要求是对称正定矩阵，且对角元素是正值，而非对角元素是非正值（即，要么是0值，要么是负值）。虽然这并不是代数多重网格的必要条件，然而AMG最初就是针对这种类型的矩阵提出来的，如果矩阵跟这种类型的矩阵偏差较远的话，那么很可能AMG并没有什么效果。

跟前面的几何多重网格算法一样，需要做下面的事情（注：虽然粗细网格并不像结构网格那样单元长度为h和2h，但是这里仍然用h和2h分别表示细和粗网格层次）

• 选取粗网格节点，一般是在原细网格节点的基础上，利用矩阵$\mathbf{A}$的系数，根据一定的条件来选取一些网格点作为粗网格；
• 插值矩阵${\mathbf{I}}_{2h}^h$，用来将粗网格上的修正值传递到细网格上；
• 限制矩阵${\mathbf{R}}_h^{2h}$，用来将细网格上的残差传递到粗网格上；
• 粗网格上的系数矩阵${\mathbf{A}}_{2h}$。

2.1 插值矩阵${\mathbf{I}}_{2h}^h$

先不管粗网格是如何得到的，假设已经整好了粗网格，那么怎么把粗网格上的值传递到细网格上呢？即，怎么由粗网格节点值插值得到细网格节点值呢？由于此时节点关系比较复杂，所以并不能一下子就把插值矩阵给整出来，而是需要逐个节点去处理才好的。

假设已经找好了粗网格节点（即从原本的细网格中筛选好了粗网格节点），即，把最初的细网格节点 { 1 , 2 , 3 , . . . , n } \{1,2,3,...,n\} {1,2,3,...,n}分为了$C\cup F$，粗网格节点集合$C$（从原本的细网格节点中筛选出来的作为粗网格的那些节点们，它们也与原来的某些细网格节点重合，也属于细网格节点），和仅仅是细网格的节点集合$F$（即从细网格中筛掉的没资格做粗网格的那些节点们，你们这些渣渣就不要掺合粗网格的计算了）。那么插值要解决的问题就是，知道了在粗网格节点$C$上的修正值$e_i$，如何将它们传递到仅仅是细网格的那些节点$F$上去呢？

如果细网格上的节点$i$是和粗网格节点重合的（$i\in C$），即之前选出来做粗网格的那些细网格节点，那么直接用$e_i$就好了，即$e_i=e_i$。

如果细网格上的节点$i$是只属于细网格的（$i\in F$），即之前选粗网格节点时被筛掉的那些细网格节点，那就需要插值来计算了，即$e_i=\sum _{j\in C_i}{\omega }_{ij}{e}_j$，注意插值的时候只能用到粗网格节点的值，而且只能用到节点$i$的邻近节点的值，还得是那些对节点$i$的影响较大的节点的值（人轻言微，影响过小的话咱们直接就把它们咔咔了，另行处理），一句话，只能用到那些对节点$i$影响较大的邻近粗网格节点的值，也就是加和范围中的集合$C_i$。

也就是说，插值矩阵所提供的插值算法应使得每个细网格上的节点$i$的修正值为
$$({\mathbf{I}}_{2h}^h\mathbf{e}{)}_i=\{\begin{array}{rcl}{e}_i& & ifi\in C\\ \sum _{j\in C_i}{\omega }_{ij}{e}_j& & ifi\in F\end{array}$$
问题来了，这个$j\in C_i$到底是哪些节点呢？${\omega }_{ij}$又该如何计算呢？

对于某个仅仅是细网格的节点$i\in F$，可以把细网格上关于它的系数矩阵$\mathbf{A}$中的那一行挑出来，并用其来反映各节点修正值的关系，即
$$a_{ii}{e}_i\approx -\sum _{j\in N_i}{a}_{ij}{e}_j$$
其中$N_i$表示在细网格中节点$i$的邻近节点，那么这些节点要么属于$C$（粗网格节点），要么属于$F$（仅为细网格节点），再考虑到节点之间影响程度的强弱，可以$N_i$分成三类来处理：

• 对节点$i$影响较强的邻近粗网格节点，即$C_i$，用于做节点$i$插值的粗网格节点集合；
• 对节点$i$影响较强的邻近细网格节点，用$D_i^s$表示（s即strong）；
• 对节点$i$影响较弱的邻近网格节点，用$D_i^w$表示（w即weak），这些节点既有可能是来自粗网格$C$的，也有可能是来自仅细网格$F$的；

那啥叫影响强，啥叫影响弱呢？可以简单地用系数矩阵$\mathbf{A}$的系数来定义两个节点的影响程度：

---

给定一个阈值因子$0<\theta \le 1$，若
− a i j ≥ θ max ⁡ k ≠ i { − a i k } -a_{ij} \ge \theta \displaystyle \max_{k\ne i}\{-a_{ik}\} −aij​≥θk​=imax​{−aik​}
则说$j$节点的值$u_j$对$i$节点的值$u_i$的影响较大，或者反过来，$i$节点的值$u_i$受（被）$j$节点的值$u_j$的影响较大，$u_i$对$u_j$的依赖较大。

---

用更通俗易懂的话来讲一下，系数$-a_{ij}$（注意$j\ne i$）（前面说过，系数矩阵的非对角元素是0或者负值，所以加了负号就是正值了）反映的是节点$j$对节点$i$的影响程度的大小，那么从第$i$行所有的系数中（除掉对角元$a_{ii}$）选择最大的那个 max ⁡ k ≠ i { − a i k } \displaystyle\max_{k\ne i}\{-a_{ik}\} k​=imax​{−aik​}，它反映的就是对节点$i$影响最强的那个节点$j$的影响程度，再乘上定义好的系数$\theta$，作为临界值 θ max ⁡ k ≠ i { − a i k } \theta \displaystyle \max_{k\ne i}\{-a_{ik}\} θk​=imax​{−aik​}，再把那些影响程度大过此临界值的节点，即 − a i j ≥ θ max ⁡ k ≠ i { − a i k } -a_{ij} \ge \theta \displaystyle \max_{k\ne i}\{-a_{ik}\} −aij​≥θk​=imax​{−aik​}的节点，选择为对$i$节点影响较大的节点，就理所当然了。显然，这个影响较强的节点跟粗网格节点没有必然联系，即，这些节点有的可能会被选成粗网格节点，而有的则不会被选成粗网格节点。

继续回来看$a_{ii}{e}_i\approx -\sum _{j\in N_i}{a}_{ij}{e}_j$，由于把$N_i$分成了三类$C_i,D_i^s,D_i^w$，所以可以把它写成
$$\begin{array}{rl}{a}_{ii}{e}_i& \approx -\sum _{j\in N_i}{a}_{ij}{e}_j\\ & =-\sum _{j\in C_i}{a}_{ij}{e}_j-\sum _{j\in D_i^s}{a}_{ij}{e}_j-\sum _{j\in D_i^w}{a}_{ij}{e}_j\end{array}$$
跟插值公式$e_i=\sum _{j\in C_i}{\omega }_{ij}{e}_j$相对照，不难发现，右端第1项就不用管了，直接把$-a_{ij}/a_{ii}$给到${\omega }_{ij}$就好了。而右端的第2和第3项是需要处理的，先看第3项，由于它们对于$i$节点的影响较弱，所以直接把它们加到对角元素上就好，即把$e_j$替换成$e_i$，或者说，把它们挪到左端去，加和到对角元上，如下：
$$\left(a_{ii}+\sum _{j\in D_i^w}{a}_{ij}\right)e_i\approx -\sum _{j\in C_i}{a}_{ij}{e}_j-\sum _{j\in D_i^s}{a}_{ij}{e}_j$$
可能有人会好奇，为啥不直接把它们给扔掉呢？因为要保持权系数加和是1的准则啊，原本的$a_{ii}=\sum _{j\in N_i}{a}_{ij}$，即$1=\sum _{j\in N_i}{a}_{ij}/a_{ii}$，如果直接把这些影响小的项抛掉，就不满足这个原则了，即得到的加权系数${\omega }_{ij}$加起来也不是$1$了。

接下来处理第2项，由于对$i$节点影响较大，这些节点的值必须考虑，但是这些节点的值又不知道，所以稍微麻烦些……把这些$D_i^s$中节点$j$的值$e_j$想办法用$C_i$中的节点值来表示，从$j$节点的离散方程中寻找答案，它们同样可以写成和$i$节点相同的形式
$$\begin{array}{rl}{a}_{jj}{e}_j& \approx -\sum _{k\in N_j}{a}_{jk}{e}_k\end{array}$$
尽管$j$节点有很多邻近节点$N_j$，但我们仅考虑哪些既属于$C_i$又属于$N_j$的那些节点，即，这些节点一方面是节点$j$的邻近节点，另一方面它们又是对节点$i$影响较大的邻近粗网格节点。用这些节点值的加权插值获取$j$节点的值，加权系数用$a_{jk}$来获取，即
$$e_j\approx \frac{\sum _{k\in C_i}{a}_{jk}{e}_k}{\sum _{k\in C_i}{a}_{jk}}$$
这样子，就可以得到粗网格向细网格插值的式子了。
$$\left(a_{ii}+\sum _{j\in D_i^w}{a}_{ij}\right)e_i\approx -\sum _{j\in C_i}{a}_{ij}{e}_j-\sum _{j\in D_i^s}{a}_{ij}\left(\frac{\sum _{k\in C_i}{a}_{jk}{e}_k}{\sum _{k\in C_i}{a}_{jk}}\right)$$
与插值式子$e_i=\sum _{j\in C_i}{\omega }_{ij}{e}_j$比照即可得到插值系数${\omega }_{ij}$的值。

光说不练假把式，需要一个小例子来搞搞明白，假设二维问题，结构网格（只是为了说明问题），均分为$n\times n$网格，每个内部节点的离散框架为
$$\left[\begin{array}{ccc}-\frac{1}{2}& -2& -\frac{1}{2}\\ -1& \frac{29}{4}& -1\\ -\frac{1}{8}& -2& -\frac{1}{8}\end{array}\right]$$
对于每个内部节点$i$，咱们用到了它的东、西、南、北、东北、西北、东南、西南共8个邻近网格点来构造离散框架。假设粗糙网格直接选择$i$节点的东、西、南、北四个节点来构造（$i$节点不是粗网格点），选择阈值系数$\theta =0.2$来定义影响较强的节点，则根据 − a i j ≥ θ max ⁡ k ≠ i { − a i k } = 0.2 ∗ 2 = 0.4 -a_{ij} \ge \theta \displaystyle \max_{k\ne i}\{-a_{ik}\}=0.2*2=0.4 −aij​≥θk​=imax​{−aik​}=0.2∗2=0.4，可知东、西、南、北这4个节点属于对$i$节点影响较强的粗网格点$C_i$，东北、西北节点属于对$i$节点影响较强的细网格节点$D_i^s$，而东南、西南节点则属于对$i$节点影响较弱的细网格节点$D_i^w$（这个例子中没有对$i$节点影响较弱的粗网格节点，不然也放在$D_i^w$中）。如图（图中对网格节点的编码是先$x$再$y$的顺序）

图中的空心圆圈标识粗网格节点，实心圆圈标识细网格节点，实线表示对$i$节点影响较强的粗网格节点$C_i$，虚线表示对$i$节点影响较强的细网格节点$D_i^s$，点线表示对$i$节点影响较弱的细网格节点$D_i^w$。

那么写出$i$节点的离散格式
$$\begin{array}{rl}\frac{29}{4}{e}_i& =\underset{C_i}{\underbrace{2e_{i+n}+2e_{i-n}+e_{i+1}+e_{i-1}}}\\ & +\underset{D_i^s}{\underbrace{\frac{1}{2}{e}_{i+n-1}+\frac{1}{2}{e}_{i+n+1}}}\\ & +\underset{D_i^w}{\underbrace{\frac{1}{8}{e}_{i-n-1}+\frac{1}{8}{e}_{i-n+1}}}\end{array}$$
$D_i^w$中的节点为东南$i-n+1$、西南$i-n-1$节点，对$i$节点影响较弱，直接将其用$e_i$代替，归并到左端项中，即
$$\begin{array}{rl}\left(\frac{29}{4}-\frac{1}{8}-\frac{1}{8}\right)e_i& =\underset{C_i}{\underbrace{2e_{i+n}+2e_{i-n}+e_{i+1}+e_{i-1}}}\\ & +\underset{D_i^s}{\underbrace{\frac{1}{2}{e}_{i+n-1}+\frac{1}{2}{e}_{i+n+1}}}\end{array}$$
$D_i^s$中的节点为东北$i+n+1$、西北$i+n-1$节点，它们对$i$节点影响较强。东北$i+n+1$节点的邻近节点中属于$i$节点的$C_i$集合的（i节点的东、西、南、北节点），只有i节点的北$i+n$、东$i+1$两节点（分别是$i+n+1$节点的西、南节点），它们在东北节点离散方程中对应的系数分别为$-1$和$-2$；而西北$i+n-1$节点的邻近节点中属于$C_i$的，只有$i$节点的北$i+n$、西$i-1$节点（$i+n-1$节点的东、南节点），它们在西北节点离散方程中对应的系数分别为$-1$和$-2$，故可将$D_i^s$中的节点，即东北$i+n+1$、西北$i+n-1$节点值表示为：
$$e_{i+n+1}=\frac{-1(e_{i+n+1}{)}_W-2(e_{i+n+1}{)}_S}{-1-2}=\frac{-1e_{i+n}-2e_{i+1}}{-3}$$
$$e_{i+n-1}=\frac{-1(e_{i+n-1}{)}_E-2(e_{i+n-1}{)}_S}{-1-2}=\frac{-1e_{i+n}-2e_{i-1}}{-3}$$
将其代入上式，可得
$$\begin{array}{rl}\left(\frac{29}{4}-\frac{1}{8}-\frac{1}{8}\right)e_i& =\underset{C_i}{\underbrace{2e_{i+n}+2e_{i-n}+e_{i+1}+e_{i-1}}}+\underset{D_i^s}{\underbrace{\frac{1}{2}{e}_{i+n-1}+\frac{1}{2}{e}_{i+n+1}}}\Rightarrow \\ 7e_i& =\underset{C_i}{\underbrace{2e_{i+n}+2e_{i-n}+e_{i+1}+e_{i-1}}}+\underset{C_i}{\underbrace{\frac{1}{2}\frac{-1e_{i+n}-2e_{i-1}}{-3}+\frac{1}{2}\frac{-1e_{i+n}-2e_{i+1}}{-3}}}\Rightarrow \\ 7e_i& =\frac{7}{3}{e}_{i+n}+2e_{i-n}+\frac{4}{3}{e}_{i+1}+\frac{4}{3}{e}_{i-1}\end{array}$$
最终得到了插值公式
$$e_i=\frac{7}{21}{e}_{i+n}+\frac{6}{21}{e}_{i-n}+\frac{4}{21}{e}_{i+1}+\frac{4}{21}{e}_{i-1}$$
可见，其满足系数加和为1的准则，且仅用到了那些对$i$节点影响较大的粗网格节点$C_i$的值，完美哈。

这时候，有人可能会想到了，对于那些$D_i^s$中的节点（对$i$节点影响较大的细网格节点），如果在它们的邻近网格节点中找不到属于$C_i$（对$i$节点影响较大的粗网格节点）的，那上面这个方法岂不是要玩完了么？没错，就是这样的，但是不会出现这个问题，因为在筛选粗网格节点的时候根据一定的准则来筛选，让选出来的粗网格节点不会出现上面这种尴尬的无从处理的情况不就好了么，下面咱们就看看怎么去筛选粗网格节点。

2.2 粗网格构造（筛选）

粗网格节点是从原本的细网格节点中筛选出来的，那么哪些节点才有资格选中作为粗节点呢？先把那些对细网格中$i$节点影响较强的节点集合定义为$S_i$（那些个对$i$节点爱的死去活来的粉丝们），再把$i$节点影响强烈的节点集合定义为$S_i^T$（那些个$i$节点爱的死去活来的明星们），注意这俩集合可未必一样啊，就像萝卜爱青菜、青菜爱大豆一样，未必是相互的。那么在粗网格节点的选择中，要满足下述两个准则H1和H2

• H1 对于细网格中的每个节点$i$，对$i$节点影响较强的节点$j\in S_i$，要么选中作为对$i$节点影响较强的粗网格节点集合$C_i$（如前所述，这个集合将用于对$i$节点的插值计算），要么（这时候它属于对$i$节点影响较强的仅细网格节点$F_i$中的$D_i^s$）其至少被一个$C_i$中的节点强烈影响着（这样子的话，就不会出现上节最后提到的那种尴尬的情况了）。
• H2 粗网格节点集$C$应该是原网格集合中满足这个特性的最大子集合，即，$C$中任何一个节点都不会强烈依赖于$C$中的其它节点。

H2还是蛮好理解的，它意思是说，选择出来的粗网格节点呢，应该是这样子的，每一个节点都不受其它节点的强烈影响，反过来说，任何一个粗网格节点对其他的粗网格节点并不会产生强烈的影响。为啥要满足这个呢？因为呢，粗网格节点之间不影响的话（影响较弱），粗网格节点就只好影响细网格节点了，这不正是我们想要的么？正好可以把粗网格上的误差修正很好地传递到细网格节点上去呢。

H1相对比较难理解一些，它的意思是说，在细网格中，某个节点$j$如果对节点$i$的影响较强，那么（1）这个节点$j$要么被选择成为粗网格节点，如果这样子的话，它就成了对$i$节点影响较强的粗网格节点集合$C_i$中的一员，在计算中做粗网格向细网格插值的时候，主要就是用到这个粗网格点上的数据了；（2）如果不幸地，这个网格节点$j$没被选择成为粗网格节点，这时候它实际上是仅细网格节点（属于集合$F$），且是对$i$节点影响较强的细网格节点（属于集合$D_i^s$），那么它至少得有一个可以依赖的属于$C_i$的节点才行，即，至少得有一个对$i$节点影响较强的粗网格节点，其同时对这个$j$节点的影响也较强，如此一来，才能在上节讲的插值过程中，把这个$j\in D_i^s$节点值表示成$k\in C_i$的形式啊，才不会出现上节课末尾所讲的那种尴尬的情况呢。

不难想象，要同时满足H1和H2是蛮难的一件事情。鉴于插值过程的思路是需要满足H1的，所以H1是无论如何都必须要严格满足的，不然没办法搞插值了，而H2则是个指导性的准则，未必需要严格满足，因为即便是多上那么几个节点，也不会对插值造成什么影响，只不过是增加了一些计算量罢了。

粗网格的筛选也并非一蹴而就的，一般要分两步来走。第一步是对网格点做个初筛（着色，coloring）的工作，即大体上把最初的细网格划分为粗网格$C$和仅细网格$F$，第二步则更为细化，将交换某些$F$和$C$中的节点，以便让H1准则严格满足。

2.2.1 第1步 初筛

一开始，给每个细网格节点$i$赋个值，用来衡量其成为粗网格节点$C$的可能性的大小，一个最简单的办法是计算那些$i$节点影响较强的节点的数目${\lambda }_i$，这些个$i$节点影响较强的节点集合是$S_i^T$（前面定义过了）。举个栗子，如果某节点对于其它7个节点影响较强，那么其成为粗网格的可能行就是7，而如果某节点对其他2个节点影响较强，其成为粗网格的可能性就是2。显然这个值越大，其越有可能成为粗网格节点集合$C$中的一员。一旦确定好了每个节点的可能性，选择具有最大可能性${\lambda }_i$的节点作为第1个粗网格节点$C$。

接下来，我们选好的这个粗节点$i$，其强烈影响着周围的节点，那么这些被其强烈影响的节点就不应该再被选入粗节点了（因为插值的时候它们可以用粗网格节点值来获取了），它们应该被放入仅细节点$F$中，即把$S_i^T$中的节点全部放入$F$中。

随后，再看看其它节点中，有没有那些对这些新的仅细节点$F$影响较强的节点，它们有可能作为新的粗网格节点呢（即$F$中的节点一方面受$i$节点的强烈影响，一方面又受到这些新节点的强烈影响，所以用它们来给$F$中的节点做插值真是再好不过了）。所以，对于每个新的仅细网格节点$F$（同时也是属于$S_i^T$，$i$节点影响较强的那些节点），即$j\in S_i^T$节点，如果有其它节点$k$强烈影响着$j$节点的话，即$k\in S_j$，把$k$节点的权值${\lambda }_k$加1，表示其成为粗网格的希望又更进了一步。简单来说，对于每个$j\in S_i^T$，将$k\in S_j$的权值${\lambda }_k$增加1。

紧接着，继续选第2个节点了，当然是选择具有最大可能性$\lambda$的节点了，作为新的粗网格点$i$放入$C$中；之后，再把这个点强烈影响的节点集合$S_i^T$中的节点作为新的唯细网格点放入$F$中；再对每个$j\in S_i^T$，考察对$j$节点影响较强的节点集合$S_j$，把每个$k\in S_j$的节点权值${\lambda }_k$加1。

如此继续进行下去，直到所有所有节点都被划分完毕，要么属于粗网格$C$，要么属于仅细网格$F$为止。

此时，需要一个例子，来使得初筛过程更加明晰。为简单起见，还是用结构网格来做示例，2维问题，离散格式为
$$\frac{1}{h^2}\left[\begin{array}{ccc}-1& -1& -1\\ -1& 8& -1\\ -1& -1& -1\end{array}\right]$$
那么对于这个例子而言，定义较强影响的阈值因子$\theta$是毫无意义的，因为周围节点的影响都是一样一样的。依据该离散格式给每个节点赋初值，反映其对其他节点的影响能力，显然，对于内部节点，其影响到周围8个节点，所以它的值是8；对于边节点，其对周围5个节点有影响，它的值是5；而对于角点，其影响着周围3个节点，它的值是3，如图，图中的连线表示节点间的相互影响，由于这个问题的离散框架是对称的，所以这些影响也是相互的，即你影响我，我也影响你，所以连线上并未用箭头来表示影响的方向性。

每个节点对其它节点的影响程度都赋值后，我们来选择影响能力最大的节点作为第一个粗网格节点，扫描的过程按照从左到右从上到下的顺序进行，那么内部节点的额影响能力最大，都是8，我们第一个选中的节点是左下方的内部节点，将其作为第1个粗网格节点，标记为红色，如图。

接下来，把第1个粗网格节点影响较强的节点$S_i^T$选作仅细网格点$F$，标记为仅细网格点，本例中为红色节点的周围8个节点，用空心圆圈标出，如图

随后，对于这些新的仅细网格点，考量每个节点又能影响到周围的哪些个新节点，并将这些新节点的值加1。换种说法是，在剩下没有分类的那些节点中，如果受到了新的唯细网格点的影响，就将它们的值加1，如图，图中的蓝色节点的值已经发生了改变。

继续上述过程，再选择第2个粗网格节点，选择剩下节点中的最大值的那个作为新的粗网格点，用红色标出；将第2个粗网格点的影响节点集合作为新的唯细网格点集合，用黑色空心圆圈标出；将剩余节点中，新的唯细网格节点的影响节点的值加1，用蓝色表示，如下图。

继续该过程，选择第3个粗网格节点，如下图，同样，红色表示第3个粗网格点，黑色圆圈表示新的唯细网格点，蓝色表示新的唯细网格点影响到的那些剩余网格的影响值发生了增加。

继续该过程，选择第4个粗网格节点，如图，红色表示第4个粗网格点，黑色圆圈为新的唯细网格点，已经没有剩余网格来更新影响值了。至此，网格初筛已经完成了。

由于这个例子比较特殊，所以在完成了第1步的网格初筛之后呢，划分的粗网格$C$和唯细网格$F$已经能够满足判据H1和H2了，所以不需要再做第2步的细化工作了。

但是，运气不会经常这么好的，有时候会出现做完了初筛工作后，判据H1不满足足的情况。咱们故意构造个例子来说明这个问题，还是上面的例子，上面的离散框架，不过把边界条件改成周期边界条件，即${\varphi }_{i,j}={\varphi }_{i\pm n,j\pm n}$，咱们这有5行5列节点，所以$n=5$，由于是周期边界条件，所以每个节点算下来影响值都是一样的，都是8，如下图。注意，图中最上方的那行和最右端的那列是为了表明是周期边界条件，其节点跟最下面的那行和最左边那列是一样的，只是为了便于理解，特此说明。

按照上述初筛算法去选择粗网格节点，那么第1个粗网格节点（左下角点），其周围唯细网格点，仅细网格点周围的影响值增加点，如下图

第2个粗网格节点（最下面一行左起第3个点），如下图。

第3个粗网格节点，（从下向上第3行，从左到右第1列的节点），如图

第4个粗网格点，（从下向上第3行，从左向右第3列的节点），如图

已经没有剩余节点了，那么把原来的细网格节点初筛成了粗网格节点$C$和仅细网格节点$F$，但是，非常明显地，这种划分并不符合准则H1，即最上面两行和最右边两列的仅细网格节点是不对的，即下图中的紫色节点是不符合H1准则的，这些点的依赖关系仅发生在唯细网格之间，无法完成后续的插值计算。

比如，以从下向上第4行从左到右第1列的节点为例（下图中黄色节点），其插值的时候需要用到周围8个节点的值（图中与黄色节点有连线的节点），在这8个节点中，只有南部S节点是粗网格点（红色的那个），而其它点全是唯细网格点，要命的是，这些7个唯细节点中，只有东南、西南、东、西这4个节点跟这个南部S节点有依赖关系，而其余的3个节点（东北、北、西北节点）全部跟这个唯一的粗网格点（南S节点）毫无依赖，这就没办法去完成插值计算了！即，初筛的结果并不满足准则1，因为准则1告诉我们，原本细网格中对节点$i$影响较强的那些节点（图中其周围的8个关联节点），要么被划分成粗网格，要么要跟被划分的粗网格产生依赖关系，而本例中的黄色节点的东北、北、西北三个节点，就属于，既没有被划分成粗网格，也没有跟被划分成粗网格的点（注意，这些点还得是$i$节点的影响节点，因为我们插值的时候只考虑$i$节点的影像节点，不再向外拓展的）发生依赖关系，所以就没办法做后续插值了。

怎么办呢？这时候就得做第2步操作，来细选网格点，让那些不符合H1准则的唯细网格点$F$点转化成粗网格点$C$点才好。

2.2.2 第2步 细化

这个第2步就是为了解决在初筛之后，出现不满足准则H1的情况，其实也可以叫做第2次筛选，名字随便起，只要知道它的作用就好了。

细化的方法不难，其思想就是测试之前选好的唯细网格$F$中的每个点$i$，看看是否能保证那些属于$D_i^s$的节点至少依赖上1个$C_i$中的节点，即，对于每个唯细网格节点$i$，对其影响较强的粗网格节点集合为$C_i$，对其影响较强的唯细网格节点集合为$D_i^s$，那么$D_i^s$中的节点是否能依赖上$C_i$中的哪怕1个节点来，这样插值的时候我们就能把$i$节点的值用$C_i$节点的值来表示了。

如果对于某个唯细网格$F$中的节点$i$，测试后发现对其影响较强的唯细网格节点$j$，不跟$C_i$节点（对$i$节点影响较大的粗网格节点）产生任何依赖关系，那么就强制把这个$j$节点转化成粗网格节点$C$，然后再对$i$节点的影响较强的唯细网格节点集合$D_i^s$进行逐个测试，看看是否满足，如果$D_i^s$中的所有节点都至少强烈依赖上一个$C_i$中的节点，那么就把$j$节点永久地转化成了粗网格节点。如果很不幸，还有某些在$D_i^s$中的节点跟$C_i$中的节点毫无瓜葛，那么只好把$i$节点自身转化成粗网格节点了，同时把$j$节点自身从粗网格节点再转化为唯细网格节点，即交换了$i$和$j$的属性。

如此反复继续下一个$F$节点的检测和属性更替，直到所有的$F$节点都检测满足为止。

仍旧用上节的最后一个不满足准则1的例子来作说明，看看如何做第2遍的筛选工作。

还是上面提到的那个黄色节点，即下图中的节点A，其周围的8个影响节点中的东北、北、西北三个节点，即图中的B、C、D节点不满足准则1，所以咱们先把北节点（B节点）转化成粗网格试试看，再次检测C、D节点发现它们和B节点是严重依赖的，所以挂靠上了一个对A节点影响较大的粗网格节点B，依旧用红色标记。

继续下一个$F$中的节点，即下图中的E节点，其周围的8个节点中，只有节点G没有挂靠上对E节点影响较大的粗网格节点，所以直接把F节点转化成粗网格节点试试看，这时候E节点的周围8个节点都满足了H1准则，OK，如图。

继续下一个$F$中的节点，下图的H点，对其影响较大的周围8个节点均满足准则1，无需处理。

继续下一个点，下图中的I点，其周围8个点中的J、K、L挂靠不上对I影响较大的粗网格点，那就把J节点转化成粗网格点，测试后发现K点倒是挂上J点了，但是L点还是没有挂上任何一个对I影响较大的粗网格点，所以，干脆把I节点转化成粗网格点，J节点转化回唯细网格点，这时候，不用再检测了，因为把I节点转化成了粗网格，将来插值的时候，I节点直接用的是粗网格的值，不用再构造插值式子了。

下一个点，M点，OK的，不需要处理。

好像扫描顺序出问题了，右边两列应该先折腾吧，哈哈，不去管它们了，反正把集合$F$中的节点挨个扫描一遍就是了。接下来扫描最上面一行的$F$集合中的节点啊。发现N节点的O节点不满足准则1，将其转化为粗网格后，满足了准则1，如图。

右上角点是满足的，不用再处理了。

接下来是右边两列节点的检测和处理，P节点，OK的。

接下来是Q节点的检测，其R节点不满足准则1，将R节点转化为粗网格后，满足准则1，如下图。

还剩3个节点了，对它们逐个检测后发现，都满足准则1，故第2遍筛选过程结束，至此我们将原本的细网格划分成了粗网格$C$和唯细网格节点$F$，且严格满足准则1，为插值计算提供了可靠的基础。顺便说一句，准则2可能未必满足，不过也无关紧要了，即便多选了几个粗网格点，顶多是增加了一些计算量罢了，只要插值方法可靠就行了。

2.3 粗网格限制矩阵${\mathbf{R}}_h^{2h}$

限制矩阵是用来把细网格上的残差传递到粗网格上去，直接用插值矩阵的转置就好
$${\mathbf{R}}_h^{2h}=({\mathbf{I}}_{2h}^h{)}^T$$

2.4 粗网格稀疏矩阵${\mathbf{A}}_{2h}$

粗网格上的系数矩阵，直接根据插值矩阵、限制矩阵来构造，和之前的几何多重网格一样的方法
$${\mathbf{A}}_{2h}={\mathbf{R}}_h^{2h}{\mathbf{A}}_h{\mathbf{I}}_{2h}^h$$

2.5 粗细网格循环算法

跟之前的几何多重网格算法完全一样，不再赘述。