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.pdfhtml

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

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

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

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

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

因为要把粗网格和细网格的信息相互传递转化,因此实际上须要3个矩阵来作这些事情(用下标 h h 来标识细网格,用下标 2 h 2h 来标识粗网格):app

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

1.1 插值矩阵 I 2 h h \bold I_{2h}^{h}

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

在这里插入图片描述
如何根据粗网格( 2 h 2h )上的变量 v \bold v 来插值算得细网格( h h )上的变量 u \bold u 呢?不难发现 u 2 , u 4 , u 6 u_2,u_4,u_6 是和 v 1 , v 2 , v 3 v_1,v_2,v_3 重合的节点上的变量值,因此直接让它们相等就行了。而对于 u 1 , u 3 , u 5 , u 7 u_1,u_3,u_5,u_7 这些点处的变量值,它们的位置是位于两个粗单元节点之间的(且是中点),因此用其两侧的粗单元节点值的算术平均来处理就行了(注意,边界节点的值是0),即 u 1 = v 1 / 2 u_1=v_1/2 u 3 = v 1 / 2 + v 2 / 2 u_3=v_1/2+v_2/2 u 5 = v 2 / 2 + v 3 / 2 u_5=v_2/2+v_3/2 u 7 = v 3 / 2 u_7=v_3/2 。写成矩阵形式,即
[ u 1 u 2 u 3 u 4 u 5 u 6 u 7 ] = 1 2 [ 1 2 1 1 2 1 1 2 1 ] [ v 1 v 2 v 3 ] \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \\ u_6 \\ u_7 \end{bmatrix} =\frac{1}{2} \begin{bmatrix} 1 & & \\ 2 & & \\ 1 & 1 & \\ & 2 & \\ & 1 & 1 \\ & & 2 \\ & & 1 \\ \end{bmatrix} \begin{bmatrix} v_1 \\ v_2 \\ v_3 \end{bmatrix}
从而得到了插值矩阵 I 2 h h \bold I_{2h}^{h}
I 2 h h = 1 2 [ 1 2 1 1 2 1 1 2 1 ] \bold I_{2h}^{h}=\frac{1}{2} \begin{bmatrix} 1 & & \\ 2 & & \\ 1 & 1 & \\ & 2 & \\ & 1 & 1 \\ & & 2 \\ & & 1 \\ \end{bmatrix} ide

1.2 限制矩阵 R h 2 h \bold R_h^{2h}

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

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

仍是结构网格,若是问题扩展为二维,即笛卡尔网格,那么这时候插值矩阵和限制矩阵是啥样的呢?测试

u i , j u_{i,j} v i , j v_{i,j} 分别表示密网格和粗网格上的变量,先看从粗网格到细网格转化变量的插值矩阵,跟一维问题相似,若是细网格节点和粗网格节点重合,那么直接让二者存储的变量相同,即
u 2 i , 2 j = v i , j u_{2i,2j}=v_{i,j}
若是细网格节点恰好落在粗网格两节点之间,就用粗网格这俩节点上变量平均来计算,即
u 2 i + 1 , 2 j = 1 2 ( v i , j + v i + 1 , j ) u 2 i 1 , 2 j = 1 2 ( v i , j + v i 1 , j )   u 2 i , 2 j + 1 = 1 2 ( v i , j + v i , j + 1 ) u 2 i , 2 j 1 = 1 2 ( v i , j + v i , j 1 ) u_{2i+1,2j}=\frac{1}{2}(v_{i,j}+v_{i+1,j}) \quad 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}) \quad u_{2i,2j-1}=\frac{1}{2}(v_{i,j}+v_{i,j-1})
若是细网格节点是位于四个粗网格节点之间的,那么用这四个粗网格节点上变量的平均来计算,即
u 2 i + 1 , 2 j + 1 = 1 4 ( v i , j + v i + 1 , j + v i , j + 1 + v i + 1 , j + 1 )   u 2 i + 1 , 2 j 1 = 1 4 ( v i , j + v i + 1 , j + v i , j 1 + v i + 1 , j 1 )   u 2 i 1 , 2 j + 1 = 1 4 ( v i , j + v i 1 , j + v i , j + 1 + v i 1 , j + 1 )   u 2 i 1 , 2 j 1 = 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}) \\ ~\\ u_{2i-1,2j-1}=\frac{1}{4}(v_{i,j}+v_{i-1,j}+v_{i,j-1}+v_{i-1,j-1}) \\
如此,即可推得插值矩阵 I 2 h h \bold I_{2h}^{h}

二维问题的限制矩阵 R h 2 h = 1 4 ( I 2 h h ) T \bold R_h^{2h}=\frac{1}{4}(\bold I_{2h}^{h})^T ,粗网格的每一个节点值由细网格上的9个节点值插值获得,各节点系数如图所示。
在这里插入图片描述


v i , j = 1 16 [ 4 u 2 i , 2 j + 2 ( u 2 i + 1 , 2 j + u 2 i 1 , 2 j + u 2 i , 2 j + 1 + u 2 i , 2 j 1 ) + ( u 2 i + 1 , 2 j + 1 + u 2 i + 1 , 2 j 1 + u 2 i 1 , 2 j + 1 + u 2 i 1 , 2 j 1 ) ] 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})]
这样,也可推得限制矩阵 R h 2 h \bold R_h^{2h}

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

1.4 粗网格系数矩阵 A 2 h \bold A_2h

粗网格上的系数矩阵 A 2 h \bold A_{2h} 是由细网格上的系数矩阵 A h \bold A_h 和变量传递关系转化而来的,故有 A 2 h = R h 2 h A h I 2 h h \bold A_{2h}=\bold R_h^{2h} \bold A_h \bold I_{2h}^{h} ,完美。

好比,若是一维问题的控制方程和零边值条件为
d 2 ϕ d x 2 = f ( x ) , 0 x 1 ϕ ( x = 0 ) = ϕ ( x = 1 ) = 0 -\frac{d^2\phi}{dx^2}=f(x),\quad 0\le x\le 1\\ \phi(x=0)=\phi(x=1)=0
跟前面讲的网格剖分方式同样,细网格等分8个单元,单元长度为 h = 1 / 8 h=1/8 ,用中心格式近似2阶项,有
d 2 ϕ d x 2 = ϕ i + 1 2 ϕ i + ϕ i 1 h 2 -\frac{d^2\phi}{dx^2}=-\frac{\phi_{i+1}-2\phi_{i}+\phi_{i-1}}{h^2}
则细网格上系数矩阵 A h \bold A_h
A h = 1 h 2 [ 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 ] \bold A_h=\frac{1}{h^2}\begin{bmatrix} 2 & -1 \\ -1 & 2 & -1 \\ & -1 & 2 & -1 \\ &&-1 & 2 & -1 \\ &&&-1 & -2 & -1 \\ &&&&-1 & 2 & -1 \\ &&&&&-1 & 2 \\ \end{bmatrix}
粗网格等分4个单元,则粗网格上系数矩阵为
A 2 h = R h 2 h A h I 2 h h = 1 4 [ 1 2 1 1 2 1 1 2 1 ] 1 h 2 [ 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 ] 1 2 [ 1 2 1 1 2 1 1 2 1 ] = 1 ( 2 h ) 2 [ 2 1 1 2 1 1 2 ] \begin{aligned} \bold A_{2h}&=\bold R_h^{2h} \bold A_h \bold I_{2h}^{h} \\ &=\frac{1}{4} \begin{bmatrix} 1 & 2 & 1 & & & & \\ & &1 & 2 & 1 & & \\ & & & &1 & 2 & 1 \end{bmatrix} \frac{1}{h^2}\begin{bmatrix} 2 & -1 \\ -1 & 2 & -1 \\ & -1 & 2 & -1 \\ &&-1 & 2 & -1 \\ &&&-1 & -2 & -1 \\ &&&&-1 & 2 & -1 \\ &&&&&-1 & 2 \\ \end{bmatrix}\frac{1}{2} \begin{bmatrix} 1 & & \\ 2 & & \\ 1 & 1 & \\ & 2 & \\ & 1 & 1 \\ & & 2 \\ & & 1 \\ \end{bmatrix} \\ &=\frac{1}{(2h)^2} \begin{bmatrix} 2& -1 \\ -1 & 2 & -1 \\ & -1 & 2 \end{bmatrix} \end{aligned}

1.5 多重网格算法流程

若是仅有两层网格,细网格-粗网格-细网格这样作完一个循环,就是两层网格的V型循环,其算法流程以下:

  1. 在细网格上迭代求解 A h u = b h \bold A_h \bold u = \bold b_h ,用Jacobi或Gauss-Seidel迭代,注意只作3次迭代步,获得还没有收敛的解 u h \bold u_h (虽然没收敛,可是高频份量衰减了不少);
  2. 计算细网格上面的残差 r h = b h A h u h \bold r_h=\bold b_h - \bold A_h \bold u_h ,并将 r h \bold r_h 传递到粗网格上面去,用 r 2 h = R h 2 h r h \bold r_{2h}=\bold R_h^{2h}\bold r_h 获得粗网格上面的残差 r 2 h \bold r_{2h}
  3. 计算粗网格上系数矩阵 A 2 h = R h 2 h A h I 2 h h \bold A_{2h}=\bold R_h^{2h} \bold A_h \bold I_{2h}^{h} (这个能够放到最前面算好存好拿来直接用,不用每次都算的),求解粗网格上方程 A 2 h E 2 h = r 2 h \bold A_{2h} \bold E_{2h}=\bold r_{2h} ,一样用迭代方法来求解,仍是迭代3步就好,迭代初值选择为 E = 0 \bold E = \bold 0 ,迭代完成后算出了修正值 E 2 h \bold E_{2h}
  4. 将粗网格上的修正值 E 2 h \bold E_{2h} 插值回细网格,用 E h = I 2 h h E 2 h \bold E_{h}=\bold I_{2h}^{h}\bold E_{2h} ,得到细网格上的修正值 E h \bold E_h
  5. 在细网格上继续求解 A h u = b h \bold A_h \bold u = \bold b_h ,仍是用迭代方法迭代3步,注意这时候迭代初值选择为 u h + E h \bold u_h+\bold E_h ,即包含了粗网格的修正成分。

多重网格通常会构造好多层的粗细网格,而粗细网格遍历循环方式能够由V、W、F型,如图
在这里插入图片描述

1.6 例子

若是没有例子,那么仍是搞不明白,这时候能够找个很简单的小例子来试试看,就用前面提到的一维零边值问题,控制方程和边界条件以下
d 2 ϕ d x 2 = f ( x ) , 0 x 1 ϕ ( x = 0 ) = ϕ ( x = 1 ) = 0 -\frac{d^2\phi}{dx^2}=f(x),\quad 0\le x\le 1\\ \phi(x=0)=\phi(x=1)=0
用前面提到的网格剖分方式,细网格等分8单元,粗网格等分4单元,插值矩阵、限制矩阵、细网格系数矩阵、粗网格系数矩阵前面都推出来了。咱们让方程中的 f ( x ) = sin ( 2 π x ) f(x)=-\sin(2\pi x) ,那么这个问题的精确解就很容易获得是 ϕ ( x ) = sin ( 2 π x ) / ( 2 π ) 2 \phi(x)=-\sin(2\pi x)/(2\pi)^2 ,即细网格上的精确解
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 \bold 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. 选择初始值 u = 0 \bold u=\bold 0 ,Gauss-Seidel迭代方法迭代3步,求解 A h u = b h \bold A_h \bold u = \bold b_h ,获得 u h \bold u_h
    u h = [ 0.0122 , 0.0172 , 0.0122 , 0.0000 , 0.0122 , 0.0172 , 0.0122 ] T \bold u_h=[ -0.0122, -0.0172, -0.0122, -0.0000, 0.0122, 0.0172, 0.0122]^T
  2. 计算细网格上面的残差 r h = b h A h u h \bold r_h=\bold b_h - \bold A_h \bold u_h ,得
    r h = [ 0.2500 , 0.3536 , 0.2500 , 0.0000 , 0.2500 , 0.3536 , 0.2500 ] T \bold r_h = [ -0.2500, -0.3536, -0.2500, -0.0000, 0.2500, 0.3536, 0.2500]^T
    r h \bold r_h 传递到粗网格上面去,用 r 2 h = R h 2 h r h \bold r_{2h}=\bold R_h^{2h}\bold r_h 获得粗网格上面的残差 r 2 h \bold r_{2h}
    r 2 h = [ 0.3018 , 0.0000 , 0.3018 ] T \bold r_{2h}=[-0.3018, -0.0000, 0.3018]^T
  3. 求解粗网格上方程 A 2 h E 2 h = r 2 h \bold A_{2h} \bold E_{2h}=\bold r_{2h} ,Gauss-Seidel迭代方法迭代3步,注意迭代初值选择为 E = 0 \bold E = \bold 0 ,迭代完成后算出了修正值 E 2 h \bold E_{2h}
    E 2 h = [ 0.0094 , 0.0000 , 0.0094 ] T \bold E_{2h}=[-0.0094, -0.0000, 0.0094]^T
  4. 将粗网格上的修正值 E 2 h \bold E_{2h} 插值回细网格,用 E h = I 2 h h E 2 h \bold E_{h}=\bold I_{2h}^{h}\bold E_{2h} ,得到细网格上的修正值 E h \bold E_h
    E h = [ 0.0047 , 0.0094 , 0.0047 , 0.0000 , 0.0047 , 0.0094 , 0.0047 ] T \bold E_h=[ -0.0047, -0.0094, -0.0047, -0.0000 , 0.0047, 0.0094, 0.0047]^T
  5. 在细网格上继续求解 A h u = b h \bold A_h \bold u = \bold b_h ,注意这时候迭代初值选择为 u h + E h \bold u_h+\bold E_h
    u h = u h + E h = [ 0.0169 , 0.0267 , 0.0169 , 0.0000 , 0.0169 , 0.0267 , 0.0169 ] T \bold u_h=\bold u_h+\bold E_h=[ -0.0169, -0.0267, -0.0169, -0.0000, 0.0169 , 0.0267 , 0.0169]^T
    显而易见,粗网格的修正起了做用,这时候的 u h \bold u_h 更加接近精确解了,采用Gauss-Seidel迭代方法再迭代3步,得
    u h = [ 0.0189 0.0257 0.0189 0.0000 0.0189 0.0257 0.0189 ] T \bold u_h=[ -0.0189, -0.0257, -0.0189, -0.0000, 0.0189, 0.0257, 0.0189]^T
    其残差 r h = b h A h u h \bold r_h=\bold b_h - \bold A_h \bold u_h 的2范数为0.2165。还能够继续进行该细-粗-细的循环直到解收敛。
    注意,若是不用多重网格,只用Gauss-Seidel迭代方法只在细网格上进行迭代,进行了6步迭代后,所获得的的值 u h \bold u_h
    u h = [ 0.0165 0.0233 0.0165 0.0000 0.0165 0.0233 0.0165 ] T \bold u_h=[ -0.0165, -0.0233 , -0.0165 , -0.0000 , 0.0165 , 0.0233 , 0.0165]^T
    其残差 r h = b h A h u h \bold r_h=\bold b_h - \bold A_h \bold u_h 的2范数为0.2500,跟多重网格相比(0.2165)仍是有差距的,代表多重网格方法的确是能起到加速收敛的做用。

2 代数多重网格(AMG)算法

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

注意,系数矩阵 A \bold A 仍然要求是对称正定矩阵,且对角元素是正值,而非对角元素是非正值(即,要么是0值,要么是负值)。虽然这并非代数多重网格的必要条件,然而AMG最初就是针对这种类型的矩阵提出来的,若是矩阵跟这种类型的矩阵误差较远的话,那么极可能AMG并无什么效果。

跟前面的几何多重网格算法同样,须要作下面的事情(注:虽然粗细网格并不像结构网格那样单元长度为h和2h,可是这里仍然用h和2h分别表示细和粗网格层次)

  • 选取粗网格节点,通常是在原细网格节点的基础上,利用矩阵 A \bold A 的系数,根据必定的条件来选取一些网格点做为粗网格;
  • 插值矩阵 I 2 h h \bold I_{2h}^h ,用来将粗网格上的修正值传递到细网格上;
  • 限制矩阵 R h 2 h \bold R_h^{2h} ,用来将细网格上的残差传递到粗网格上;
  • 粗网格上的系数矩阵 A 2 h \bold A_{2h}

2.1 插值矩阵 I 2 h h \bold I_{2h}^h

先无论粗网格是如何获得的,假设已经整好了粗网格,那么怎么把粗网格上的值传递到细网格上呢?即,怎么由粗网格节点值插值获得细网格节点值呢?因为此时节点关系比较复杂,因此并不能一会儿就把插值矩阵给整出来,而是须要逐个节点去处理才好的。

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

若是细网格上的节点 i i 是和粗网格节点重合的( i C i \in C ),即以前选出来作粗网格的那些细网格节点,那么直接用 e i e_i 就行了,即 e i = e i e_i=e_i

若是细网格上的节点 i i 是只属于细网格的( i F i \in F ),即以前选粗网格节点时被筛掉的那些细网格节点,那就须要插值来计算了,即 e i = j C i ω i j e j e_i=\displaystyle \sum_{j \in C_i}\omega_{ij}e_j ,注意插值的时候只能用到粗网格节点的值,并且只能用到节点 i i 的邻近节点的值,还得是那些对节点 i i 的影响较大的节点的值(人轻言微,影响太小的话我们直接就把它们咔咔了,另行处理),一句话,只能用到那些对节点 i i 影响较大的邻近粗网格节点的值,也就是加和范围中的集合 C i C_i

也就是说,插值矩阵所提供的插值算法应使得每一个细网格上的节点 i i 的修正值为
( I 2 h h e ) i = { e i i f   i C j C i ω i j e j i f   i F (\bold I_{2h}^h \bold e)_i = \left\{ \begin{array}{rcl} e_i & & if~i \in C\\ \displaystyle \sum_{j \in C_i}\omega_{ij}e_j&& if ~ i \in F \end{array} \right.
问题来了,这个 j C i j \in C_i 究竟是哪些节点呢? ω i j \omega_{ij} 又该如何计算呢?

对于某个仅仅是细网格的节点 i F i \in F ,能够把细网格上关于它的系数矩阵 A \bold A 中的那一行挑出来,并用其来反映各节点修正值的关系,即
a i i e i j N i a i j e j a_{ii} e_i \approx - \sum_{j \in N_i} a_{ij} e_j
其中 N i N_i 表示在细网格中节点 i i 的邻近节点,那么这些节点要么属于 C C (粗网格节点),要么属于 F F (仅为细网格节点),再考虑到节点之间影响程度的强弱,能够 N i N_i 分红三类来处理:

  • 对节点 i i 影响较强的邻近粗网格节点,即 C i C_i ,用于作节点 i i 插值的粗网格节点集合;
  • 对节点 i i 影响较强的邻近细网格节点,用 D i s D_i^s 表示(s即strong);
  • 对节点 i i 影响较弱的邻近网格节点,用 D i w D_i^w 表示(w即weak),这些节点既有多是来自粗网格 C C 的,也有多是来自仅细网格 F F 的;

那啥叫影响强,啥叫影响弱呢?能够简单地用系数矩阵 A \bold A 的系数来定义两个节点的影响程度:


给定一个阈值因子 0 < θ 1 0\lt\theta\le1 ,若
a i j θ max k i { a i k } -a_{ij} \ge \theta \displaystyle \max_{k\ne i}\{-a_{ik}\}
则说 j j 节点的值 u j u_j i i 节点的值 u i u_i 的影响较大,或者反过来, i i 节点的值 u i u_i 受(被) j j 节点的值 u j u_j 的影响较大, u i u_i u j u_j 的依赖较大。


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

继续回来看 a i i e i j N i a i j e j a_{ii} e_i \approx - \displaystyle\sum_{j \in N_i} a_{ij} e_j ,因为把 N i N_i 分红了三类 C i , D i s , D i w C_i,D_i^s,D_i^w ,因此能够把它写成
a i i e i j N i a i j e j = j C i a i j e j j D i s a i j e j j D i w a i j e j \begin{aligned} 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{aligned}
跟插值公式 e i = j C i ω i j e j e_i=\displaystyle \sum_{j \in C_i}\omega_{ij}e_j 相对照,不难发现,右端第1项就不用管了,直接把 a i j / a i i -a_{ij}/a_{ii} 给到 ω i j \omega_{ij} 就行了。而右端的第2和第3项是须要处理的,先看第3项,因为它们对于 i i 节点的影响较弱,因此直接把它们加到对角元素上就好,即把 e j e_j 替换成 e i e_i ,或者说,把它们挪到左端去,加和到对角元上,以下:
( a i i + j D i w a i j ) e i j C i a i j e j j D i s a i j e j \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 i i = j N i a i j a_{ii}=\displaystyle \sum_{j \in N_i} a_{ij} ,即 1 = j N i a i j / a i i 1=\displaystyle \sum_{j \in N_i} a_{ij}/a_{ii} ,若是直接把这些影响小的项抛掉,就不知足这个原则了,即获得的加权系数 ω i j \omega_{ij} 加起来也不是 1 1 了。

接下来处理第2项,因为对 i i 节点影响较大,这些节点的值必须考虑,可是这些节点的值又不知道,因此稍微麻烦些……把这些 D i s D_i^s 中节点 j j 的值 e j e_j 想办法用 C i C_i 中的节点值来表示,从 j j 节点的离散方程中寻找答案,它们一样能够写成和 i i 节点相同的形式
a j j e j k N j a j k e k \begin{aligned} a_{jj} e_j & \approx - \sum_{k \in N_j} a_{jk} e_k \end{aligned}
尽管 j j 节点有不少邻近节点 N j N_j ,但咱们仅考虑哪些既属于 C i C_i 又属于 N j N_j 的那些节点,即,这些节点一方面是节点 j j 的邻近节点,另外一方面它们又是对节点 i i 影响较大的邻近粗网格节点。用这些节点值的加权插值获取 j j 节点的值,加权系数用 a j k a_{jk} 来获取,即
e j k C i a j k e k k C i a j k e_j \approx \frac{\displaystyle\sum_{k\in C_i}a_{jk}e_k}{\displaystyle\sum_{k\in C_i}a_{jk}}
这样子,就能够获得粗网格向细网格插值的式子了。
( a i i + j D i w a i j ) e i j C i a i j e j j D i s a i j ( k C i a j k e k k C i a j k ) \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{\displaystyle\sum_{k\in C_i}a_{jk}e_k}{\displaystyle\sum_{k\in C_i}a_{jk}} \right)
与插值式子 e i = j C i ω i j e j e_i=\displaystyle \sum_{j \in C_i}\omega_{ij}e_j 比照便可获得插值系数 ω i j \omega_{ij} 的值。

光说不练假把式,须要一个小例子来搞搞明白,假设二维问题,结构网格(只是为了说明问题),均分为 n × n n\times n 网格,每一个内部节点的离散框架为
[ 1 2 2 1 2 1 29 4 1 1 8 2 1 8 ] \begin{bmatrix} -\displaystyle\frac{1}{2} & -2 & -\displaystyle\frac{1}{2} \\ -1 & \displaystyle\frac{29}{4} & -1 \\ -\displaystyle\frac{1}{8} & -2 & -\displaystyle\frac{1}{8} \end{bmatrix}
对于每一个内部节点 i i ,我们用到了它的东、西、南、北、东北、西北、东南、西南共8个邻近网格点来构造离散框架。假设粗糙网格直接选择 i i 节点的东、西、南、北四个节点来构造( i i 节点不是粗网格点),选择阈值系数 θ = 0.2 \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 ,可知东、西、南、北这4个节点属于对 i i 节点影响较强的粗网格点 C i C_i ,东北、西北节点属于对 i i 节点影响较强的细网格节点 D i s D_i^s ,而东南、西南节点则属于对 i i 节点影响较弱的细网格节点 D i w D_i^w (这个例子中没有对 i i 节点影响较弱的粗网格节点,否则也放在 D i w D_i^w 中)。如图(图中对网格节点的编码是先 x x y y 的顺序)

在这里插入图片描述

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

那么写出 i i 节点的离散格式
29 4 e i = 2 e i + n + 2 e i n + e i + 1 + e i 1 C i + 1 2 e i + n 1 + 1 2 e i + n + 1 D i s + 1 8 e i n 1 + 1 8 e i n + 1 D i w \begin{aligned} \frac{29}{4}e_i&=\underbrace{2e_{i+n}+2e_{i-n}+e_{i+1}+e_{i-1}}_{C_i}\\ &+\underbrace{\frac{1}{2}e_{i+n-1}+\frac{1}{2}e_{i+n+1}}_{D_i^s}\\ &+\underbrace{\frac{1}{8}e_{i-n-1}+\frac{1}{8}e_{i-n+1}}_{D_i^w} \end{aligned}
D i w D_i^w 中的节点为东南 i n + 1 i-n+1 、西南 i n 1 i-n-1 节点,对 i i 节点影响较弱,直接将其用 e i e_i 代替,归并到左端项中,即
( 29 4 1 8 1 8 ) e i = 2 e i + n + 2 e i n + e i + 1 + e i 1 C i + 1 2 e i + n 1 + 1 2 e i + n + 1 D i s \begin{aligned} \left( \frac{29}{4}-\frac{1}{8}-\frac{1}{8} \right) e_i&=\underbrace{2e_{i+n}+2e_{i-n}+e_{i+1}+e_{i-1}}_{C_i}\\ &+\underbrace{\frac{1}{2}e_{i+n-1}+\frac{1}{2}e_{i+n+1}}_{D_i^s} \end{aligned}
D i s D_i^s 中的节点为东北 i + n + 1 i+n+1 、西北 i + n 1 i+n-1 节点,它们对 i i 节点影响较强。东北 i + n + 1 i+n+1 节点的邻近节点中属于 i i 节点的 C i C_i 集合的(i节点的东、西、南、北节点),只有i节点的北 i + n i+n 、东 i + 1 i+1 两节点(分别是 i + n + 1 i+n+1 节点的西、南节点),它们在东北节点离散方程中对应的系数分别为 1 -1 2 -2 ;而西北 i + n 1 i+n-1 节点的邻近节点中属于 C i C_i 的,只有 i i 节点的北 i + n i+n 、西 i 1 i-1 节点( i + n 1 i+n-1 节点的东、南节点),它们在西北节点离散方程中对应的系数分别为 1 -1 2 -2 ,故可将 D i s D_i^s 中的节点,即东北 i + n + 1 i+n+1 、西北 i + n 1 i+n-1 节点值表示为:
e i + n + 1 = 1 ( e i + n + 1 ) W 2 ( e i + n + 1 ) S 1 2 = 1 e i + n 2 e i + 1 3 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 = 1 ( e i + n 1 ) E 2 ( e i + n 1 ) S 1 2 = 1 e i + n 2 e 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}
将其代入上式,可得
( 29 4 1 8 1 8 ) e i = 2 e i + n + 2 e i n + e i + 1 + e i 1 C i + 1 2 e i + n 1 + 1 2 e i + n + 1 D i s 7 e i = 2 e i + n + 2 e i n + e i + 1 + e i 1 C i + 1 2 1 e i + n 2 e i 1 3 + 1 2 1 e i + n 2 e i + 1 3 C i 7 e i = 7 3 e i + n + 2 e i n + 4 3 e i + 1 + 4 3 e i 1 \begin{aligned} \left( \frac{29}{4}-\frac{1}{8}-\frac{1}{8} \right) e_i&=\underbrace{2e_{i+n}+2e_{i-n}+e_{i+1}+e_{i-1}}_{C_i}+\underbrace{\frac{1}{2}e_{i+n-1}+\frac{1}{2}e_{i+n+1}}_{D_i^s} \Rightarrow \\ 7 e_i&=\underbrace{2e_{i+n}+2e_{i-n}+e_{i+1}+e_{i-1}}_{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}}_{C_i} \Rightarrow \\ 7 e_i&=\frac{7}{3}e_{i+n}+2e_{i-n}+\frac{4}{3}e_{i+1}+\frac{4}{3}e_{i-1} \end{aligned}
最终获得了插值公式
e i = 7 21 e i + n + 6 21 e i n + 4 21 e i + 1 + 4 21 e i 1 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 i 节点影响较大的粗网格节点 C i C_i 的值,完美哈。

这时候,有人可能会想到了,对于那些 D i s D_i^s 中的节点(对 i i 节点影响较大的细网格节点),若是在它们的邻近网格节点中找不到属于 C i C_i (对 i i 节点影响较大的粗网格节点)的,那上面这个方法岂不是要玩完了么?没错,就是这样的,可是不会出现这个问题,由于在筛选粗网格节点的时候根据必定的准则来筛选,让选出来的粗网格节点不会出现上面这种尴尬的无从处理的状况不就行了么,下面我们就看看怎么去筛选粗网格节点。

2.2 粗网格构造(筛选)

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

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

H2仍是蛮好理解的,它意思是说,选择出来的粗网格节点呢,应该是这样子的,每个节点都不受其它节点的强烈影响,反过来讲,任何一个粗网格节点对其余的粗网格节点并不会产生强烈的影响。为啥要知足这个呢?由于呢,粗网格节点之间不影响的话(影响较弱),粗网格节点就只好影响细网格节点了,这不正是咱们想要的么?正好能够把粗网格上的偏差修正很好地传递到细网格节点上去呢。

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

不难想象,要同时知足H1和H2是蛮难的一件事情。鉴于插值过程的思路是须要知足H1的,因此H1是不管如何都必需要严格知足的,否则没办法搞插值了,而H2则是个指导性的准则,未必须要严格知足,由于即使是多上那么几个节点,也不会对插值形成什么影响,只不过是增长了一些计算量罢了。

粗网格的筛选也并不是一蹴而就的,通常要分两步来走。第一步是对网格点作个初筛(着色,coloring)的工做,即大致上把最初的细网格划分为粗网格 C C 和仅细网格 F F ,第二步则更为细化,将交换某些 F F C C 中的节点,以便让H1准则严格知足。

2.2.1 第1步 初筛

一开始,给每一个细网格节点 i i 赋个值,用来衡量其成为粗网格节点 C C 的可能性的大小,一个最简单的办法是计算那些 i i 节点影响较强的节点的数目 λ i \lambda_i ,这些个 i i 节点影响较强的节点集合是 S i T S_i^T (前面定义过了)。举个栗子,若是某节点对于其它7个节点影响较强,那么其成为粗网格的可能行就是7,而若是某节点对其余2个节点影响较强,其成为粗网格的可能性就是2。显然这个值越大,其越有可能成为粗网格节点集合 C C 中的一员。一旦肯定好了每一个节点的可能性,选择具备最大可能性 λ i \lambda_i 的节点做为第1个粗网格节点 C C

接下来,咱们选好的这个粗节点 i i ,其强烈影响着周围的节点,那么这些被其强烈影响的节点就不该该再被选入粗节点了(由于插值的时候它们能够用粗网格节点值来获取了),它们应该被放入仅细节点 F F 中,即把 S i T S_i^T 中的节点所有放入 F F 中。

随后,再看看其它节点中,有没有那些对这些新的仅细节点 F F 影响较强的节点,它们有可能做为新的粗网格节点呢(即 F F 中的节点一方面受 i i 节点的强烈影响,一方面又受到这些新节点的强烈影响,因此用它们来给 F F 中的节点作插值真是再好不过了)。因此,对于每一个新的仅细网格节点 F F (同时也是属于 S i T S_i^T i i 节点影响较强的那些节点),即 j S i T j \in S_i^T 节点,若是有其它节点 k k 强烈影响着 j j 节点的话,即 k S j k \in S_j ,把 k k 节点的权值 λ k \lambda_k 加1,表示其成为粗网格的但愿又更进了一步。简单来讲,对于每一个 j S i T j \in S_i^T ,将 k S j k \in S_j 的权值 λ k \lambda_k 增长1。

紧接着,继续选第2个节点了,固然是选择具备最大可能性 λ \lambda 的节点了,做为新的粗网格点 i i 放入 C C 中;以后,再把这个点强烈影响的节点集合 S i T S_i^T

相关文章
相关标签/搜索