- A translation to: Chapter 7-3 of Strang, Gilbert. Computational Science and Engineering
- Translator: Adversarr
考虑如下的两点边值问题:
$$
\begin{cases}
- u''(x) = f(x)\quad \forall x \in (0, 1) \
u(0) = u(1) = 0
\end{cases}
$$
我们使用有限差分法进行求解,对于空间采用$N+1$点均匀划分:
$$
x_0 = 0, ..., x_n = i / N, ..., x_N = 1
$$
二阶导数用如下的公式近似:
$$
f(x_j) = u''(x_j) \approx \frac{2u(x_{j}) - u(x_{j-1}) - u(x_{j+1})}{(1/N)^2}
$$
这转化为求解如下的线性方程组求解问题:
$$
\frac{1}{h^2} \begin{bmatrix}
2 & -1 & 0 & 0& \cdots & 0\
-1 & 2 & -1& 0& \cdots & 0\
0 & -1 & 2 & -1& \cdots& 0\
0 & 0& -1 & 2& \cdots& 0\
\vdots& \vdots&\vdots&\vdots & \ddots & \vdots\
0 & 0 & 0 & 0&\cdots & 2
\end{bmatrix} \begin{bmatrix}
u(x_1)\u(x_2)\u(x_3)\
\vdots\
u(x_{N-1})
\end{bmatrix} = \begin{bmatrix}
f(x_1)\f(x_2)\f(x_3)\
\vdots\
f(x_{N-1})
\end{bmatrix}
$$
我们将方程组的求解问题写成
Note:这里使用有限差分法是可以收敛到精确解的,可以考虑在使用等距网格的情况下的有限元方法,它们诱导了相同的线性方程组求解问题。
我们一般将这里的求解选取的
Jacobi 和 Gauss- Seidel 迭代都能够产生“光滑”的误差:误差向量
Background:考虑线性方程组
$A x = b$ 的求解,使用迭代求解。
$x_0$ 为迭代初值;$e_i = A^{-1}b - x_i$ ,表示当前解与真解的误差;$r_i = Ae_i = b - A x_i$ 表示当前解表达方程组的误差Multi Grid方法最常见的可以用于求解Poisson问题的解。均匀网格下的有限差分法将这样的Poisson问题变成线性方程组的求解问题。
多重网格法(MultiGrid)可以极大的解决迭代求解器的效率问题,其思想是将原本网格上的线性方程求解问题转化到一个更粗的网格上,在更粗的网格上施用迭代求解器(例如Jacobi/GS求解器),细网格上的低频分量转化为了粗网格上的高频分量,粗网格求解器消除的高频分量实则就是细网格上的低频分量。
Recall:例如间隔为$\Delta x$的网格上$e^{i\omega \cdot n \Delta x}$ 分量,在$2 \Delta x$ 上对应$e^{i 2 \omega \cdot n \Delta x}$,频率也提升了一倍。
在粗网格上,原本的低频分量很快就被“干掉”了。而且,我们在粗网格上只需要很少次数的迭代,就能够使得这些分量消失。这样的性质使得多重网格法能够使用固定次数的迭代,就能较精确地求解现实中很多的稀疏方程组(这里的次数不随着矩阵规模而变大)。
多重网格法一般用于对称的矩阵求解问题。它的最核心要素是用两个矩阵
-
限制矩阵(Restriction Matrix)
$R=R_{2h}^h$ :把向量从细网格传到粗网格上; - 插值矩阵(Interpolation matrix)$I= I_{2h}^h$:将粗网格上的结果传到细网格上
- 原始矩阵在粗网格上的近似
$A_{2h}=RA_hI$
Note:这里用
$h, 2h$ 这样的记号是参考了有限差分求解的网格粗细记号,用来表征从一个粒度为$h$ 的细网格,和另一个粒度为$2h$ 的粗网格。
我们先考虑这些矩阵的定义:假设粗网格上有 3 个点,
$$
v = (v_1, v_2, v_3)^T
$$
我们希望得到细网格上的数值,可以用线性插值得到,写成矩阵形式可以得到上面定义的插值矩阵如下。注意到$u$ 的2/4/6行恰好是
Note:利用两点边值问题的边界条件。
考虑
Note:这一段说的有点奇怪,简而言之就是,在细网格上的任何一个低频误差分量都能在粗网格上有一一的体现。(Recall:高频误差分量都能够被迭代求解器很好的处理)
反过来,我们可以从细网格限制到粗网格上(类似于一个“平均”)。它将$u$从一个细网格传递到粗网格上。一种可能的是直接令
二维空间中,粗网格到细网格的插值是通过双线性插值得到。双线性插值可以通过$x, y$两个方向上分别进行线性插值得到。
水平方向的插值:
$u_{2i, 2j} = v_{i, j}$ $u_{2i, 2j + 1} = \frac 1 2 (v_{i, j} + v_{i, j+1})$
竖直方向的插值:
$u_{2i +1, 2j} = (v_{i, j} + v_{i+1, j}) / 2$ $u_{2i +1, 2j+1} = (v_{i, j} + v_{i+1, j} + v_{i, j+1} + v_{i+1, j+1}) / 4$
我们不给出插值算子
Recall:有限维线性空间到有限维线性空间的线性映射,选取一个合适的基的情况下可以表达为一个矩阵。
反过来,从细网格到粗网格的算子
其实你不难发现这些稀疏的由来,原本一维情况下加权平均系数是$(\frac14, \frac12, \frac14)$。$R2D$在这一个局部上的系数正是这个小的向量的外积。而$R2D$正是$R$与自身的Kronecker Product。
Note:考虑 3x3 的网格和,标号为
1 2 3 4 5 6 7 8 9
它拍平后为
1 2 3 4 5 6 7 8 9
Kronecker积恰好对应了这样的变换。简单验证下,$R$为3x7的矩阵,$R2D$ 是 9x49 的矩阵,维数也是相对应的。
Recall:矩阵的 Kronecker 积对应了 Tensor Product
现在,我们能够在网格之间相互传递向量,这就是几何理解的多重网格法,我们能在网格基于间距h,2h和4h的情况下进行计算。这个想法可以延伸到三角形元素(每个三角形自然地分成四个三角形)。实际计算的几何形状可以比我们在一个正方形网格考虑的模型更加复杂。
当几何上的多重网格变得较为复杂,或者我们只是获得了一个矩阵,而没有相应的“网格”概念,这是我们可以使用Algebraic multigrid(AMG)。它也是基于多尺度求解的思想,可以直接应用在$Ax = b$问题上,而不依赖于具体的网格。
Jacobi Iteration: $$ x_{n+1} = D^{-1}(b - R x_n)\quad \text{where}\quad A = D + R $$
我们先考虑只使用两个网格进行求解的过程。在每个网格内的迭代,我们考虑使用 Jacobi 格式或者GS格式求解。我们前面多次指出,对于一个较大规模的问题,在细网格上的迭代法对于低频误差去除得很慢(也就是原文说的“对于精确解的低频分量收敛很慢”)。多重网格法将细网格的残差$r_h = b - A u_h$传到粗网格上$r_{2h}$;然后在粗$2h$网格上求解$A_{2h} E_{2h} = r_{2h}$;然后将$E_{2h}$传到细网格上为
这样的循环就是一个两网格的 V-cycle。我们将其称之为一个 v-cycle(小v),以下是它的步骤(残差是指
-
Iterate(迭代):在细网格上迭代,求解
$A_h u_h = b_h$ (例如使用3次的Jacobi或GS迭代) - Restrict(限制):将残差传递到粗网格上,$r_{2h} = R_h^{2h} r_h$
-
Solve(求解):在粗网格上迭代,求解修正项
$A_{2h} E_{2h} = r_{2h}$ (例如,从$E=0$开始3次迭代,在更进一步的方法中,这个方程的求解方法是任意的) -
Interpolate(插值):修正项从粗网格上传到细网格上,$E_h = I_{2h}^h E_{2h}$,将
$E_h$ 加到$u_h$ 上 -
Iterate(迭代):以
$u_h + E_h$ (消除了部分低频分量)为初值迭代,类似第一步的过程
Note:
- 残差$r$是指模型观测值与真实值的误差,模型观测值指的是
$Ax$ ,真实值就是$b$ ,$r = b - Ax$- 误差$e$是指当前解和真实解的误差,$e = A^{-1} b - x$
- 误差和残差相比,要求知道真实解的值(知道真解的值时,已经解决了矩阵求解问题)但残差只需要计算$b-Ax$,这里面每一项都是已知的。
我们希望修正项
$E_h = I_{2h}^h E_{2h}$ 能够尽可能好的表达$e_h$ ,即$E_h \approx e_h$ ,这样:$$ \begin{aligned} &u_h + I_{2h}^h E_{2h} \approx u_h + e_h = A^{-1}b\ \implies & A_hI_{2h}^{h} E_{2h}\approx b - A_hu_h = r_h\ \implies & \underbrace{R_{h}^{2h} A_h I_{2h}^h}{A{2h}} E_{2h} \approx R_{h}^{2h} r_h \end{aligned} $$
上面的最后一步给出了多重网格第三步的解释。
2-3-4 步给出的 限制-粗解-插值 过程是多重网格法的核心。回忆我们定义的三个矩阵: $$ \begin{aligned} A &= A_{h} = \text{original matrix}\ R &= R_{h}^{2h} = \text{restriction matrix}\ I &= I_{2h}^h = \text{interpolation matrix}\ \end{aligned} $$
第三步需要第四个矩阵$A_{2h}$,可以从上面的式子中定义(如下)。
例如,细网格上有7个数据点,粗网格有3个数据点,我们求得的矩阵恰好是$3\times 3$的。
在一维情况下,$A = A_h$ 可能是一个二阶差分矩阵
对于$A_{2h}$的一个自然的选择是(直接用前文中的插值矩阵$I_{2h}^h$ 代入也是相同的结果): $$ A_{2h} = RAI = \frac{1}{(2h)^2} \begin{bmatrix} 2 & -1\ -1& 2 \end{bmatrix} $$ 不难看出,这恰好是$2h$的网格上的二阶差分矩阵(因为网格为$2h$宽,分母自动放大)。
Note:差分算子K$的定义从两点边值问题的矩阵求解形式中得到。
Recall:之前定义的两点边值问题,Poisson方程(微分方程)被转化为线性方程组(差分方程)。写出网格$h = 1/3$的矩阵,恰好是这里的$A_{2h}$。
Note:我的理解是,用相同的数值方法,但用不同的网格密度求解同一个方程,粗网格的解$u_{2h}$插值到细网格上,应该“接近”细网格的解$u_h$。
我们选取的$R = \frac{1}{2} I^T$还有一个很好的性质:如果原先的矩阵
读者自证不难。
Note:对于Poisson 问题的有限元法,
- 除去纯Neumann边值问题,其他的情况下由双线性形式诱导的矩阵都是正定的;
- 纯Neumann边值问题添加$\int_{\Omega} u \mathrm d x = 0$ 的条件后,诱导的矩阵也是正定的;
尽管第一步和第五部不是多重网格法的核心,但第一步和第五步是必要的:第一步(smoother)将原本误差中的高频分量消除,而第五步(post-smoother)将粗网格插值引入的一小部分高频误差消除。在实现上,各种Jacobi和GS求解器都是可行的。
Note:这里杂糅了一些自己的理解,原文为
Steps 1 and 5 are necessary, but they are really outside the essential multigrid idea. The smoother is step 1, the post-smoother is step 5. Those are normal iterations for which weighted Jacobi or Gauss-Seidel is satisfactory.
假设在第三步我们精确地求出了粗网格上的误差方程。多重网格误差修正项$E_h$是否能等于
Note:看之前的步骤部分后面给的一段分析。
在这里,方便起见,在不引起歧义的情况下省略了一部分上下标(公式打得我好累)
作用在$e$上的四个矩阵:
- 第二步,残差计算涉及了$e$:$r_h = b- A_h u_h = Ae$;
- 第二步,传到粗网格时,限制矩阵作用在其上:$r_{2h} = RA_he$
- 第三步,假设精确地求解$2h$上的方程:$E_{2h} = A_{2h}^{-1}r_{2h} = (RAI)^{-1}R A_h e$
- 第四步,传回细网格,插值矩阵作用在其上,最终形式如下
Note:不能将$A_{2h}^{-1}=(RAI)^{-1}$拆开,这应该很显然。
当
Note:列空间上恒等的说明。$S$的列空间是说 ${ Su: u \in \mathbb R^n },在列空间中的任一向量
$Su$ ,因为$S$幂等,所以 $$ S(Su) = S^2 u = Su $$ 这就说明$S$在列空间上的矩阵是恒等矩阵。
Note:其实有一个挺常见的结论,矩阵若满足幂等条件
$S^2 = S$ ,那么$S$是投影矩阵。(并不是正交投影)
因此,多重网格法求出的的修正项
Note: 第一步和第五步分析实际上是指对于Jacobi、GS格式的分析,迭代法的误差通常都能写成 $$ e_{n+1} = e_n - P^{-1} A e_n $$ 的形式,例如 Jacobi 迭代($D = \mathrm{diag} A$): $$ x_{n+1} = D^{-1}(b - R x_n)\quad \text{where}\quad A = D + R $$ “简单”计算可得: $$ e_{n+1} = (I - D^{-1} A) e_n $$
我们直接给出之前例子的$S$矩阵如下:
该矩阵的特征值是
Note:这里的“这样的向量不会出现在粗网格上”原文为 “Those are vectors that don’t appear on the coarse grid.”
作者表达的意思并不是说这些误差向量不会体现到粗网格上,而是说粗网格是无法消除这些误差向量。这听起来很奇怪是吧?
我们换一种说法来表述这个过程。因为
Note:特征向量组就是原本空间的一组基,误差可以在这组基下分解,然后考虑各个分量在多重网格法下如何变化。
对于
对于
放下原文体会一下: The other
$5-2$ grid points account for the nullspace of$S$ , where$E = Se = 0$ means no improvement from multigrid.
备注:满足$Su = 0$的“高频分量”$(u_1, 0, u_3, 0, u_5)$并不是三个离散正弦序列(这在后面的Fourier分析会考虑)。严谨的说法是,$I_{2h}^h$ 的列空间等于
Recall:多重网格法将低频分量转化为粗网格上的高频分量,再迭代求解器消除。
如果我们的迭代求解器能够完美地消除误差的高频分量,而多重网格能够消除剩下的低频分量,那么一个v-cycle就会求出精确解!但这再现实中显然是不可能的。幸运的是如果我们仔细分析可以发现,一次v-cycle能够将误差固定地缩小一个常数倍,并且这个常数是与网格密度$h$无关的:
通常这里的
收敛因子:迭代矩阵的谱半径(a spectral radius of the overall iteration matrix)
Note:不管是插值、限制,还是Jacobi迭代、GS的迭代的一步,都是$O(n)$的,这里的$n$说的是稀疏系统的规模,例如非零元素个数。
Note:网格点数为$n$的情况下,矩阵规模是$n \times n$的(有$n^2$个自由度,这里的矩阵维数说的是$n^2$),但实际上非零元素的个数是$O(n)$的(不管是有限元还是有限差分法都是)
Recall:Jacobi迭代对于$n\times n$矩阵是需要$O(n^2)$次收敛。
关于迭代次数和精度的关系还有一点可以补充。用户可能希望数值解的误差和方程离散引入的误差相同。之前的例子中,我们的二阶差分是2阶精度的,因此希望$e = O(h^2) = O(N^{-2})$。这意味着我们多次(而非固定次数)的v-cycle来达到精度要求。为了达到
除了多次执行v-cycle,或者将他们嵌套在一起称为 V-Cycle 或 W-Cycle,一种更好的方案是使用全多重网格(Full Multigrid, i.e. FMG)。这些循环都在下一章中介绍。对于我们需要的$O(h^2)$精度,FMG方法的复杂度是$O(n)$的。
我们之前对于第三步的求解方法不加以限制。显而易见的是我们可以在原本的v-cycle第三步也执行类似的多重网格法。举个栗子,我们在第三步仅仅使用了Jacobi迭代,原本的最低频分量的频率(尽管频率提升了一倍)还是低的。但如果不断地进行v-cycle,那么这些分量能在$4h$、$8h$等等,以至于非常粗的网格上得到消除(例如$512h$)。
我们有一个很自然的选择,将多层网格划分为多层的(例如3层,左下图),它先不断加粗网格到$(2h, 4h, 8h)$上,然后直接逐层返回到$(4h, 2h, h)$上。这个嵌套的 v-cycle 序列就是一个V-cycle(大写V)。需要注意的是,在更粗的网格上的矩阵求解是更快的。一些分析指出,我们可以花费稍多一点的时间在粗网格上是值得的。我们构造了如下中图的W-cycle,它能够消除更多的低频分量,通常相较于V-cycle有更好的结果。
FMG方法“渐进”地优于的V/W-cycle,如右上图所示。它从最粗的网格出发,先算出$8h$上的解, 然后插值到$4h$上,来获得$4h$上的一个求解初值,这个初值应当接近于
Note:联系微分方程的数值解,粗网格上的解为细网格上的求解提供了良好的初值。
Note:从FMG计算的流程不难推断,其为$h$网格提供的初值,相应的低频误差很“少”。
尽管V-cycle(或其他的两个算法,但我们只分析V-cycle)总操作数多于v-cycle,但它们只相差一个常数倍。这是因为每当我们切换到更粗一倍的网格的时候,我们求解的规模也随之降低一半。对于一个维度为$d$的微分方程,粗网格上的问题规模是细网格的$1/2^d$。下面的公式大致估计了V-cycle的复杂度: $$ V < (1 + \frac{1}{2^d} + \left(\frac{1}{2^d}\right)^2 + \cdots )v < \frac{2^d}{2^d - 1} v $$ 其中的$V$表示V-cycle的操作数,$v$表示v-cycle的操作数。使用类似的记号,用$F$表示FMG的操作数,类似的有: $$ F < \frac{2^d}{2^d - 1} V < \left(\frac{2^d}{2^d - 1}\right)^2 v $$ 顺便一提的是,FMG方法在实际应用的多,但需要你有比较好的编程能力。
原本一次v-cycle的我们给出了矩阵$S$来刻画误差的变化,对于一个三层的V-cycle,我们也希望有一个相应的矩阵$S_3$来刻画误差的变化。在这里分析时,我们不考虑平滑(Smoothing)的影响。$S$和$S_3$只是将原本的问题投影到更粗的网格上而已。就其本身而言,如果没有平滑步的情况下$S_3$并不是一个好的求解器。
Note:解释一下这里的平滑(Smoothing)的含义。考虑一个三层的V-Cycle:
$h$ -Iterate$h\rightarrow 2h$ -Restrict:$2h$ -Solve(also a v-cycle):
$2h$ -Iterate Smoothing$2h\rightarrow 4h$ -Restrict$4h$ Solve$4h \rightarrow 2h$ -Interpolate$2h$ -Iterate Smoothing$2h\rightarrow h$ -Interpolate$h$ -Iterate我们在分析的时候,为了方便起见,认为误差并没有经过$2h$-Iterate的两步,否则整个分析过程会变的相当复杂。与此同时,如果我们考虑的是那些非常低频的分量(因为只有这些分量会被更粗的网格解决),$2h$-Iterate也几乎不对这些分量产生影响。
在最外层的第三步求解的问题实际上是$A_{2h}E_{2h}=r_{2h}$,我们只需要重新利用$A_{2h}^h= RA_hI$这样的公式就能得到:
其中的$R_{2h}^{4h}$将$2h$网格点限制到$4h$上,反过来,$I_{4h}^{2h}$将$4h$网格点插值到$2h$上。如果我们选择的网格密度为$h=1/16$,那么这个矩阵乘法是
和之前相同的分析,在一次不考虑平滑过程的V-cycle下,误差为$(I-S_3) e$。在$h=1/16$的网格上,一共有15个不同的频率分量,只有最低的3个被$S_3 e$(大致上)消除。我们现在仅仅是求解了一个$3\times3$的问题,考虑到在$h$网格和$2h$网格上的迭代步(Iterate),高频的和中间频率的误差分量都被这些迭代求解器消除了。
TODO: 这里省去了原文的“数值试验”这一小节。
读者可能已经意识到,一个矩阵(例如v-cycle的$I-S$)就足以描述多重网格法中的一个步骤。而对于FMG方法,尽管我们希望能得到它的特征值和特征向量,但它们通常实在是难以求出。先前的数值试验的结果让我们对多重网格法充满信心。计算和实验给我们一种分析、改进(Diagnostic)的手段,来发现在哪里收敛得慢,需要作出一些改变(通常在边界条件上需要特殊处理)。我们也需要一个先验的预测工具,最好的就是模态分析(Modal Analysis)。
本文分析的核心是观察 Fourier 谐波在多重网格求解时的变化。在我们的例子中,考虑到零编制条件,我们选取离散的正弦波来分析。我们一直使用Poisson方程导出的问题(Model Problem),来观察多重网格法矩阵$I-S$对这些正弦波进行了如何的作用。这一小节的最终结果(从数学上)说明了为何多重网格法是有效的,也说明了两种频率分量是如何被混合的。即$I-S$的特征向量是由两个谐波混合而来。
Note: 频率为$k$的正弦波是指: $$ y(i) = \sin 2 \pi \frac{ik}{N} $$ 它能满足我们的零边值条件
我们直接指出,混合在一起的分量是频率为$k$和$N-k$的(我们只考虑$k \le N-k$的情况,这样$k$始终是频率更低的)。写出两个频率的正弦波
其中
对于我们的问题(二阶差分算子矩阵求解),我们可以观察它们乘以$S = I(RAI)^{-1}RA$后的情况。但事实上,多重网格法对应的矩阵实际上是$(I-S)$,它给出了经过234步后误差和修正项之间的差$e - E = (I-S)e$。所以我们直接观察$I-S$的在低频正弦波$y_k$和高频正弦波$Y_k$上的作用:
对于v-cycle,我们能直接分析得到这个美妙的结果(Beautiful formula),但在更复杂的问题中是很难直接分析得到的。在此我们指出多重网格法的四个关键点:
- 误差的低频分量被多重网格很好地去除(例子中$k=1,2,3$):在上面的公式中,$\cos$项是很接近$1$的(因为低频分量$k\ll N$),它大致是$O((k/N)^2)$的。这说明多重网格法几乎消除了误差中所有的低频分量。这些分量就是在细网格上Jacobi/GS迭代器难以去除的部分。
- 一对正弦波被多重网格法混叠:混叠是由$R, I$矩阵带来的。它们不像$A$,$A$是方阵(Square Matrix),它能够将不同的频率分量分开(因为正弦波就是二阶差分算子的特征向量,就如正弦波是二阶导数算子的特征值一样)。对于一般的矩阵,会将不同的分量混叠在一起。
- 两个频率混合得到的$e = y_k + Y_k$是$I-S$的特征向量,对应的特征值为$1$:我们将上面的式子直接相加就能得到这个结果。因为$Y_k$和$y_k$的唯一区别就是符号,那么$y_k + Y_k$有$(e_1, 0, e_3, 0 ,...)$这样的形式。它们正是我们在之前找到的$S$的零空间中的向量。因为$Se = 0$,这些分量不被多重网格法解决。
-
$I-S$ 的其他的特征向量是$S y_k$,对应特征值为$0$:之前我们指出$S^2 = S$,因此有$(I-S) Sy_k = 0$。在我们的例子中,我们显式地求出了$Sy_1$和$Sy_2$,它们分别是$(1, 2, 2, 2, 1)$和$(1, 2, 0, -2, 1)$的倍数。它们大多是一些低频分量,并且多重网格法将这些分量从误差中去除。
根据这个公式,我们还能得到$S y_k$不过是$y_k$和$Y_k$的组合。为了找到更好的组合$e ^ $(这个$e^
这个混合频率分量$e^*$被多重网格法完全消除了(这对于所有的$k \le N-k$都成立)。
总的来说,由第一步计算后得到的光滑的误差向量$e_h$在$y_k + Y_k$这个方向上只有一些很小的分量(因为细网格上的Jacobi迭代能够有效得消除误差的高频分量)。多重网格法并不会改变这个(小的)分量。更大的分量是沿着$e^*$的,这些分量能被多重网格法完全消除。
Note:为何在$y_k + Y_k$上只有很小的分量,而沿着$e^{*}$由更大的分量?
需要注意$e^{}$的表达式中,$y_k$和$Y_k$的系数。当$k \ll N$时,有$\cos \approx 1$,那么这个$e^ \approx 2y_k$。如果考虑误差向量$y=\lambda_1 y_k + \lambda_2 Y_k$用$y_k+Y_k, e^$表示,那么 $y \approx \frac12 (\lambda_1 -\lambda_2) e^{} + \lambda_2 (y_k+Y_k)$,于此同时,我们知道误差中的高频分量少,即$\lambda_2$是很小的,这就能够印证上面的说法。
为了让以上的分析更加完整,我们下面来看看一对频率分量是如何被混叠在一起的。这样的混叠现象产生于限制矩阵
Note:这里的上标表示的是不同网格下的正弦波。
上面的公式直接指出了$R$对不同频率的混叠。我们无法仅仅从输出中分析原本输入的向量是$y_k$还是$Y_k$(方程组欠定)。另一方面,如果我们直接观察$R$的行、列数,这样的“欠定现象”对于一个矮、胖的矩阵是再常见不过的。在粗网格上的不同频率分量数目大致上只有细网格的一半。
Note:这里说的不同频率分量的数目也就是其线性空间的维数。
反过来,$R^T$对输入的向量有与$R$相反的作用。当$R$将两个分量混合为一个的时候,插值矩阵以(混合后的)这一个频率分量作为输入,输出两个频率分量:
从这个公式中我们还能够看出,对一个粗网格上的一个向量插值会额外产生一个细网格上的高频分量。但是这个高频分量的振幅很小,因为$\cos k\pi / N\approx 1$。
Recall:这说明为何算法的第五步是必要的,为消除低频误差而采用的多重网格法(第2-4步)还会额外引入一个很小的高频分量,第五步就能够将这些分量再去除掉。下面对例子的分析也说明了这一点。
我们直接指出,如果我们没有第1、5步的光滑过程,而反复执行第2-4步是没有任何作用的!考虑$e = y_k + Y_k$误差,$(I-S)e = e$ 说明了它在整个2-4步的迭代后没有产生任何变化。只有在1、5步的光滑作用(Smoothing Matrix,即前文提到的$M = I - P^{-1}A$)能够很好的将这些高频误差分量消除。
假设在第1、5步采用(单步的)加$w = \frac{2}{3}$权的Jacobi迭代,即$M = I - \frac{1}{3} A$,那么整个1-5步对应的矩阵为
Note:矩阵特征值绝对值的最大值决定了其最多能将向量放大多少倍,也能说明至少缩小多少倍。
Recall: 矩阵的谱半径
Note:$(I-S)$的秩仅仅为3,所以剩下的特征值为0,整个矩阵的谱半径就是$1/9$。
Note: 这里的矩阵$M(I-S)M$是作用在误差上的,即经过完整的1-5步后,误差$e$降低为$M(I-S)Me$。
这说明,原本$I-S$矩阵$\lambda = 1$的三个特征方向对应了$M(I-S)M$矩阵的$\lambda = 1/9$的特征方向。
如果我们不用多重网格法,$M$的最大特征值(谱半径)为$0.91$。这还只是在第1、5步只进行了一次光滑的效果(只进行了一次Jacobi迭代),采用更多次的光滑会带来更好的效果。
这一节的分析,我们忽略零边界条件,而考虑更一般的情况。与此同时,我们考虑一个无穷大的网格。之前的问题中,零边界条件让我们只考虑了正弦波,而当边界条件消失,研究对象从原本的有限长的正弦波$y_k$向量,变为无限长、任意频率的简谐波向量$y_\omega$。通常用连续频率$\omega$的复指数函数来表示(连续频率是指$\omega$不必再是整的或是有理的):
我们需要一个“无穷维的矩阵”$K_\infty$来作用在这些无穷维的向量上。二阶差分格式
对于$K_\infty$的作用的分析为我们提供了$A_h$的作用的信息。同时,当我们切换到粗网格上,$K_\infty$也告诉我们$A_{2h}$的信息。
(在这种更广义的情况下)限制矩阵$R$也会产生频率混叠,不同的是,现在混合的频率分量是$\Omega$和$\Omega + \pi$。这两个频率分量之间仅仅相差了一个符号,即$e^{i\pi n} = (-1)^n$。这与我们之前对$y_k, Y_k$的分析是一致的。
本节采用的Fourier分析自始至终都会对无穷序列$(I-S)y_\omega$和$(I-S)y_{\omega + \pi}$进行分析,整个过程在有穷维的情况是相仿的。对于有穷维情况的分析说明了为何多重网格法是有效的,这些推导过程对于无穷维的情况也是成立的。
你可能觉得这一节是不重要的,因为实际求解过程中不可能出现无穷维、没有边界的情况。但为何纯模态分析依然是重要的?这是因为这样的分析中,我们包含了所有的Fourier序列(原文:It allows Fourier to work freely)!Poisson方程$\mathrm{div}(c \cdot \nabla u) = f$(例如守恒型方程的微分形式)在更一般的区域上是相当复杂的。但如果我们将$c$限制维一个常数并且忽略所有的边界条件,这些“内部的特征向量”(Interior Eigenvectors)就是这些$y_\omega, y_{\omega + \pi}$的组合。Fourier方法并不擅长求解带边值条件的问题,这为我们也留下了一个难题,即如何正确地处理边界。
作为本章节(也是本文)的结尾,我们花一小部分介绍下代数多重网格法(Algebraic Multigrid, AMG)。先前我们从一个实际的Poisson问题出发,研究了如何对于二阶差分方程采用多重网格进行求解。更一般的$Au =b$问题中并不涉及“网格”的概念(There is no grid in the background)。原本多重网格法的一些重要概念需要在一般的问题中被定义出:光滑向量、节点的邻接关系、粗网格。前两个概念是显而易见的,但是由于我们没有网格,从$A_h$得到$A_{2h}$的方法是比较模糊的。
Note:看作图的邻接矩阵,每一个非零的$i, j$元素一定程度上表明节点$i$和$j$的“亲密程度”或是“流量”。
-
光滑向量(Smooth vectors):是指那些$| u |$和$| Au|$是同数量级(comparable)的向量。$u$中的高频分量会被$A$放大,就如同二阶微分/差分算子对于$\sin k t$放大了$k^2$倍那样。
-
节点的邻接关系(Connected nodes):在网格上,节点的邻接关系是由$A$的非零元描述的。如果没有网格,我们直接观察矩阵,那些有较大影响的(significant)元素告诉我们
$i$ 和$j$ 之间是“邻接”的(例如$|A_{ij}|> |A_{ii}|/10$ ) -
粗网格(Coarse subset of nodes
$C$ ):每一个有较大影响的$A_ij$都说明了$u_j$对于$u_i$有较大的影响。这意味着,$e_j$和$e_i$有很大概率是可比的的(comparable,同数量级的含义)。因此,对于$e_i$和$e_j$而言,我们不需要它们都出现在粗网格上。假设在细网格上两个点是邻接的,那么它们很可能不会都在出现在粗网格$C$中。例如,如果节点$i$不在粗网格$C$中,那么每一个、或至少一个对$i$有很强影响的节点$j$都应当在粗网格中。有一系列的启发式的规则、算法来选择粗网格中包含的节点。但总体上说,使用包含较多节点的粗网格通常优于包含较少节点的粗网格。
书()中还给出了AMG方法的插值函数。它从在$C$中的节点$j$的误差$E_j$开始分析,并使它们保持不变。如果$i$不在
对$e_j$的插值组合方法会将最大的权重赋予对$i$有最大影响的$j$。但较小的那些元素也不能被完全忽略。换言之,最终采用插值的方法并不是简单的平均。这也说明AMG方法在计算成本比一般的多重网格法上是更加昂贵,但它能够适合更一般的稀疏矩阵求解(并且有软件可以控制AMG的行为,不需要我们烦心)。我们仍然希望“光滑+多重网格”的强强联手能够让AMG方法在$O(n)$内收敛到$Au = b$的精确解。
在此说一些题外话,对于固体力学和流体力学的计算之间有一个显著的差异:对于固体,原本的有限元网格相对而言是较粗的,通常我们寻求的是“工程上的精确”(Engineering Accuracy),并没有很细节上的运动现象需要考虑。对结构问题而言,多重网格法并不是一个常用的方法(消元法更简单)。但流体力学不一样,流体通常有很多细小的运动(例如高雷诺数情况下的湍流)。不仅是从数值方法上还是物理含义上,多重网格法都成为了一个自然的选择。尽管由于对流项的存在、有限元方法的“迎风修正”(Upwind adjustment),一般的流体并不产生一个对称的矩阵。特别是当我们应用多尺度方法的时候,多重网格法的收敛分析随之变得更为复杂。