DOI:10.12677/AAM.2021.1012442,PDF,HTML,XML,被引量下载: 489浏览: 706国家自然科学基金支持
作者:张 衡：无损检测技术福建省高等学校重点实验室(福建技术师范学院)，福建 福清；福建技术师范学院大数据与人工智能学院，福建 福清；湖南三一工业职业技术学院工程机械学院，湖南 长沙 ；张 宏：无损检测技术福建省高等学校重点实验室(福建技术师范学院)，福建 福清；福建技术师范学院电子与机械工程学院，福建 福清；游文杰：福建技术师范学院大数据与人工智能学院，福建 福清；郑汉垣：龙岩学院传播与设计学院，福建 龙岩
关键词:病因抑制方法；病态结构；病态因子；去病因子；最优预条件子；Method of Inhibiting Pathogeny；Ill-Condition Structure；Ill-Condition Factor；Ill Elimination Factor；Optimal Preconditioner

文章引用：张衡, 张宏, 游文杰, 郑汉垣. 二维边值问题有限元离散系统病态问题的病因抑制方法[J]. 应用数学进展, 2021, 10(12): 4162-4171. https://doi.org/10.12677/AAM.2021.1012442

1. 引言

偏微分方程的数值求解问题，通常转化为离散系统——如 $Ax=b$的线性方程组求解 [1] ，离散系统的病态(高条件数， $\text{Cond}\left(A\right)$表示 $A$的条件数)问题随着求解规模增加而更严重 [2] ，成为制约大规模求解效率和精度的瓶颈因素，因此，在求解过程中，使用适当的预条件子进行预处理，将离散系统化成条件数低、易求解的同解方程，是成功求解的关键 [3] [4] 。

目前关于病态问题以及预处理技术的理论和方法仍然不完善，缺乏对病因的研究，预条件子的设计、预处理效果的评价多用数值模拟的手段进行探索和验证 [5] - [11] ，产生的预处理方法使用范围小，可移植性差，难以找到通用的、最优的预条件子 [12] ，很多学者甚至认为不存在通用的预处理方法 [13] [14] [15] 。目前关于病态问题的病因研究鲜见有成果发表。

基于病因的结构分析方法 [16] [17] [18] [19] [20] ，本文提出离散系统病态问题的病因抑制方法：确定病因，精准抑制病因发作。

针对矩形区域上变分形式的二维泊松方程边值问题，使用有限元方法求解，分别基于三角剖分和四角剖分，使用Lagrange形函数产生的离散系统，是稀疏病态方程组。本文使用病因抑制方法解决该方程的病态问题，通过结构分析 [16] [17] [18] [19] [20] ，确定病态结构，分离出病态因子、构造去病因子。对于原发病态和继发病态，分别提出评估条件数的方法和预处理策略，定量分析以去病因子作为预条件子的预处理效果(有效性、低成本)，说明了在几乎不显著增加计算量的情况下，预处理后的条件数关于离散系统规模一致有界，预处理后离散系统保持正定对称。

2. 病态问题的病因抑制方法

2.1. 病态因子与去病因子

定义1：如果矩阵 $Z\in R^{\alpha \times \beta },\beta >\alpha$， $\alpha ,\beta$是整数， $Z{Z}^{\text{T}}$可逆， $\text{Cond}\left(Z{Z}^{\text{T}}\right)$是阶数 $\alpha$的增函数。则称 $Z$为 $\alpha \times \beta$病态因子。

定义2：对于定义1中的 $\alpha \times \beta$病态因子 $Z$，如果有可逆矩阵 $H$满足： $Z{Z}^{\text{T}}=H{H}^{\text{T}}$，则称矩阵 $H$为属于 $Z$的去病因子 [16] 。

命题1 [16] ：设 $\beta \ge \alpha$， $Z$是 $\alpha \times \beta$病态因子， $P$是 $\beta$阶正定对称矩阵，则有

1). $Cond\left(ZP{Z}^{\text{T}}\right)\le Cond\left(P\right)Cond\left(Z{Z}^{\text{T}}\right)$，

2). 如果矩阵 $H$为属于 $Z$的去病因子，则有 $Cond\left(H^{-1}ZP{Z}^{\text{T}}{H}^{-\text{T}}\right)\le Cond\left(P\right)$。

2.2. 矩阵的病态结构

定义3：可逆矩阵 $A$的如下结构称为病态结构：

$A=ZP{Z}^{\text{T}}+Q$(1)

其中 $A,Q\in R^{\alpha \times \alpha }$， $P\in R^{\beta \times \beta }$， $\beta >\alpha$， $\alpha ,\beta$是整数， $Z$是 $\alpha \times \beta$病态因子， $P$是 $\beta$阶正定对称矩阵， $\Vert ZP{Z}^{\text{T}}\Vert \gg \Vert Q\Vert$， $\text{Cond}\left(P\right)$关于 $A$的阶数一致有界， $\text{Cond}\left(ZP{Z}^{\text{T}}\right)$是 $A$的阶数 $\alpha$的增函数。称 $ZP{Z}^{\text{T}}$为 $A$的病态主体；也称 $Z$是 $ZP{Z}^{\text{T}}$或者 $A$的病态因子。

在定义3中，使用 $ZP{Z}^{\text{T}}$代替 $A$，相对误差为 $\Vert A-ZP{Z}^{\text{T}}\Vert /\Vert ZP{Z}^{\text{T}}\Vert =\Vert Q\Vert /\Vert ZP{Z}^{\text{T}}\Vert \ll 1$，所以不妨认为 $A\approx ZP{Z}^{\text{T}}$， $Cond\left(A\right)\approx Cond\left(ZP{Z}^{\text{T}}\right)$， $Q$对病态的影响很小，根据命题1结论1)， $Cond\left(A\right)\approx Cond\left(ZP{Z}^{\text{T}}\right)\le Cond\left(P\right)Cond\left(Z{Z}^{\text{T}}\right)$，即可以利用病态因子的条件数，估计 $A$的病态。

定义4：在定义3中，由病态因子表达的病态称为原发性病态，其他因素表达的病态称为继发性病态。

2.3. 病态问题的病因抑制方法

对于方程组

$Ax=b$(2)

如果 $A$有(1)式的病态结构，使用属于 $Z$的去病因子 $H$作为预条件子，则方程(2)化成

$H^{-1}A{H}^{-\text{T}}{x}^{\prime }=\left(H^{-1}ZP{Z}^{\text{T}}{H}^{-\text{T}}+H^{-1}Q{H}^{-\text{T}}\right)x^{\prime }=b^{\prime }$(3)

其中 $x^{\prime }=H^{\text{T}}x$， $b^{\prime }=H^{-1}b$。

根据命题1结论2)，在式(3)中， $Cond\left(H^{-1}A{H}^{-\text{T}}\right)\approx Cond\left(H^{-1}ZP{Z}^{\text{T}}{H}^{-\text{T}}\right)\le Cond\left(P\right)$这说明属于 $Z$的去病因子 $H$抑制或者消除了病态因子 $Z$的作用，即消除了原发性病态。因为 $\text{Cond}\left(P\right)$关于 $A$的阶数一致有界，因此， $H{H}^{\text{T}}$是最优预条件子 [4] 。

综上所述，为解决方程(2)的病态问题，提出如下病因抑制方法：通过对病态矩阵 $A$进行结构分析，分离出病态因子 $Z$，确定如(1)式的病态结构；一方面可以用病态因子 $Z$转置积的条件数评价 $A$的病态程度；另一方面可以针对病态因子 $Z$(无须针对 $A$)设计最优预条件子——去病因子 $H$，用于消除 $A$的原发性病态。

病因抑制方法的关键是确定病态因子。

2.4. 一个特别的病态因子及其去病因子

设 $diag\left(S_1,\cdots ,S_n\right)$表示对角元素为矩阵 $S_i\left(1\le i\le n\right)$的块对角矩阵，记 $0_{\alpha }$是 $\alpha$维零向量， ${\omega }_{i,j}^{\left(n\right)}=\sqrt{2/\left(n+1\right)}\sin \left[ij\pi /\left(n+1\right)\right]={\omega }_{j,i}^{\left(n\right)}$， ${\Omega }_n=\left({\omega }_{i,j}^{\left(n\right)}\right)$是n阶正弦变换矩阵， ${\lambda }_{i,j}^{\left(n,m\right)}=4\sqrt{1-{\cos }^2\left[i\pi /\left(2n+2\right)\right]{\cos }^2\left[j\pi /\left(2n+2\right)\right]}$， ${\Lambda }_i=diag\left({\lambda }_{i,1}^{\left(n,m\right)},\cdots ,{\lambda }_{i,m}^{\left(n,m\right)}\right)$，

$H=\left({\Omega }_n\otimes {\Omega }_m\right)\Lambda$， $\Lambda =diag\left({\Lambda }_1,\cdots ,{\Lambda }_n\right)$

$D^{\text{T}}=\left({\sigma }_0^{\text{T}},{\sigma }_1^{\text{T}},{\sigma }_2^{\text{T}},{\sigma }_3^{\text{T}}\right)$， ${\sigma }_0=\left(1,1,-1\right)$， ${\sigma }_1=\left(-1,1,1\right)$， ${\sigma }_2=\left(1,-1,1\right)$， ${\sigma }_3=\left(-1,-1,-1\right)$

$Z=\left(I_n,0_n\right)\otimes Z_1+\left(0_n,I_n\right)\otimes Z_0$， $Z_i=\left(I_m,0_m\right)\otimes {\sigma }_{2i+1}+\left(0_m,I_m\right)\otimes {\sigma }_{2i}$， $i=0,1$，则有

命题2：1). $Z{Z}^{\text{T}}=H{H}^{\text{T}}$；2). $Cond\left(Z{Z}^{\text{T}}\right)={\left({\lambda }_{n,m}^{\left(n,m\right)}/{\lambda }_{1,1}^{\left(n,m\right)}\right)}^2$；3). $Cond\left(Z{Z}^{\text{T}}\right)=O\left(\min \left(n^2,m^2\right)\right)$。

证明：1). 记 $I_n$是n阶单位矩阵， $T_n={\left(1,2,1\right)}_n$是n阶三对角矩阵， ${\bar{\Lambda }}_n=ding\left({\bar{\lambda }}_1^{\left(n\right)},\cdots ,{\bar{\lambda }}_n^{\left(n\right)}\right)$， ${\bar{\lambda }}_i^{\left(n\right)}=2\cos i\pi /\left(2n+2\right)$，则容易验证： $Z{Z}^{\text{T}}=16I_n\otimes I_m-T_n\otimes T_m$， $T_n={\Omega }_n{\bar{\Lambda }}_n^2{\Omega }_n$[21] ，因此

$16I_n\otimes I_m-{\bar{\Lambda }}_n^2\otimes {\bar{\Lambda }}_m^2={\Lambda }_{}^2$， $T_n\otimes T_m=\left({\Omega }_n\otimes {\Omega }_m\right)\left({\bar{\Lambda }}_n^2\otimes {\bar{\Lambda }}_m^2\right)\left({\Omega }_n\otimes {\Omega }_m\right)$

$H{H}^{\text{T}}=\left({\Omega }_n\otimes {\Omega }_m\right)\left(16I_n\otimes I_m-{\bar{\Lambda }}_n^2\otimes {\bar{\Lambda }}_m^2\right)\left({\Omega }_n\otimes {\Omega }_m\right)=Z{Z}^{\text{T}}$。证毕。

2). 根据命题2结论1)， $Z{Z}^{\text{T}}$的特征值为 ${\left({\lambda }_{i,j}^{\left(n,m\right)}\right)}^2$， $1\le i\le n$， $1\le j\le m$。所以， $Cond\left(Z{Z}^{\text{T}}\right)={\lambda }_{\max }\left(Z{Z}^{\text{T}}\right)/{\lambda }_{\min }\left(Z{Z}^{\text{T}}\right)={\left({\lambda }_{n,m}^{\left(n,m\right)}/{\lambda }_{1,1}^{\left(n,m\right)}\right)}^2$。证毕。

3). 记 $g\left(x,y\right)=\left(x+y+xy\right)/\left(1+x+y\right)$， $x>1$， $y>1$， ${\alpha }_n={\cot }^2\pi /\left(2n+2\right)$，则容易验证： ${\alpha }_n=O\left(n^2\right)$， $2\min \left(x,y\right)+2\ge 2g\left(x,y\right)\ge \min \left(x,y\right)$， $Cond\left(Z{Z}^{\text{T}}\right)=g\left({\alpha }_n,{\alpha }_m\right)$，所以， $g\left(x,y\right)=O\left(\min \left(x,y\right)\right)$。结论成立。证毕。

命题3：对于任意 $3\left(m+1\right)\left(n+1\right)$阶正定对称矩阵 $P$，有 $Cond\left(ZP{Z}^{\text{T}}\right)\le Cond\left(P\right)Cond\left(Z{Z}^{\text{T}}\right)$。

证明：根据命题1结论1)可得。证毕。

根据定义1、定义2和命题2， $Z$是 $mn\times 3\left(m+1\right)\left(n+1\right)$病态因子；矩阵 $H$是属于 $Z$的去病因子。

3. 二维泊松方程边值问题的Lagrange有限元离散系统

3.1. Lagrange有限元离散格式

使用有限元方法，求解如下变分形式的二维泊松方程边值问题 [22]

${\displaystyle {\iint }_D\left[p\left(\nabla u\cdot \nabla v\right)+quv\right]dxdy}+{\displaystyle {\int }_{\partial D}\sigma puvds}-{\displaystyle {\iint }_{D}fvdxdy}=0,\forall v$(4)

其中 $\nabla u=\left(\partial u/\partial x,\partial u/\partial y\right)$， $p=p\left(x,y\right)\ge p_{\min }>0$， $q=q\left(x,y\right)\ge 0$， $\sigma =\sigma \left(x,y\right)\ge 0$， $f=f\left(x,y\right)$， $\left(x,y\right)\in G=\left[a_1,b_1\right]\times \left[a_2,b_2\right]\subset R^2$。

对G做剖分： $a_1=x_0<x_1<\cdots <x_n<x_{n+1}=b_1$， $a_2=y_0<y_1<\cdots <y_m<y_{m+1}=b_2$，记 $h_{i,x}=x_i-x_{i-1}$， $h_{j,y}=y_j-y_{j-1}$， $e_{i,j}=\left[x_{i-1},x_i\right]\times \left[y_{j-1},y_j\right]$， $u_{i,j}=u\left(x_i,y_j\right)$， $v_{i,j}=v\left(x_i,y_j\right)$。 $0\le i\le n+1$， $0\le j\le m+1$, $h_{\min }=\underset{i,j}{\min }\left\{h_{i,x},h_{j,y}\right\}$， $h=\underset{i,j}{\max }\left\{h_{i,x},h_{j,y}\right\}$，则(4)式变成

${\displaystyle \underset{j=1}{\overset{n+1}{\sum }}{\displaystyle \underset{i=1}{\overset{m+1}{\sum }}\left\{{\displaystyle {\iint }_{e_{i,j}}\left[p\left(\nabla u\cdot \nabla v\right)+quv\right]dxdy}+{\displaystyle {\int }_{\partial e_{i,j}}\sigma puvds}-{\displaystyle {\iint }_{e_{i,j}}fvdxdy}\right\}}}=0,\forall v$(5)

取 $\left(\xi ,\eta \right)$为自然坐标 [23] ， $\left(\xi ,\eta \right)\in e=\left[0,1\right]\times \left[0,1\right]$，在 $e_{ij}$上，做变换 $x=x_i\left(\xi \right)=x_{i-1}+\xi h_{i,x}$， $y=y_j\left(\eta \right)=y_{j-1}+\eta h_{j,y}$及其逆变换 $\xi =\left(x-x_{i-1}\right)/h_{i,x}$， $\eta =\left(y-y_{j-1}\right)/h_{j,y}$，记 $X_{ij}=\left(x_i\left(\xi \right),y_j\left(\eta \right)\right)$， $p_{ij}=p\left(X_{ij}\right)$， $q_{ij}=q\left(X_{ij}\right)$， ${\sigma }_{ij}=\sigma \left(X_{ij}\right)$， $f_{ij}=f\left(X_{ij}\right)$， $J_{ij}=h_{i,x}{h}_{j,y}$，则直接验证可以得到

命题4：1) $\partial \left(x_i\left(\xi \right),y_j\left(\eta \right)\right)/\partial \left(\xi ,\eta \right)=diag\left(h_{i,x},h_{j,y}\right)$；

2) $\partial \left({\xi }_i\left(x\right),{\eta }_j\left(y\right)\right)/\partial \left(x,y\right)=diag\left(h_{i,x}^{-1},h_{j,y}^{-1}\right)$；

3) $J_{ij}=\left|\partial \left(x_i\left(\xi \right),y_j\left(\eta \right)\right)/\partial \left(\xi ,\eta \right)\right|=O\left(h^2\right)$；

证明：直接验证。证毕

3.1.1. 基于四角剖分的Lagrange有限元离散格式

在 $e_{ij}$上，取2次Lagrange形函数为： $N_{\alpha ,\beta }=N_{\alpha ,\beta }\left(\xi ,\eta \right)\in P_2\left[\xi ,\eta \right]$， $\left(\xi ,\eta \right)\in e$， $\left\{\alpha ,\beta \right\}\subset \left\{0,1\right\}$，满足： $N_{\alpha ,\beta }\left(\alpha ,\beta \right)=1$， $N_{\alpha ,\beta }\left(\xi ,\eta \right)=0$， $\left(\xi ,\eta \right)\ne \left(\alpha ,\beta \right)$， $\left\{\xi ,\eta \right\}\subset \left\{0,1\right\}$。

记 ${\bar{u}}_{ij}={\left(u_{i-1,j-1},u_{i-1,j},u_{i,j-1},u_{i,j}\right)}^{\text{T}}$， ${\bar{v}}_{ij}={\left(v_{i-1,j-1},v_{i-1,j},v_{i,j-1},v_{i,j}\right)}^{\text{T}}$， $U_{ij}={\bar{u}}_{ij}^{\text{T}}N$， $V_{ij}={\bar{v}}_{ij}^{\text{T}}N$

$e={\left(0,0,0,1\right)}^{\text{T}}$， $N={\left(N_{0,0},N_{0,1},N_{1,0},N_{1,1}\right)}^{\text{T}}$， $\bar{N}=diag\left(N_{0,0},N_{0,1},N_{1,0},N_{1,1}\right)$， $M_{ij}={\nabla \tilde{N}|}_{e_{ij}}$

$\stackrel{⌢}{N}=\tilde{N}{\tilde{N}}^{\text{T}}-\frac{1}{4}\left(\begin{array}{ccc}{N}_{0,0}+N_{0,1}& N_{0,0}& N_{0,1}\\ N_{0,0}& N_{0,0}+N_{1,0}& N_{1,0}\\ N_{0,1}& N_{1,0}& N_{1,0}+N_{0,1}\end{array}\right)$， $\tilde{N}=\frac{1}{2}\left(\begin{array}{c}{N}_{0,0}+N_{0,1}\\ N_{0,0}+N_{1,0}\\ N_{1,0}+N_{0,1}\end{array}\right)=\left(\begin{array}{c}{\tilde{N}}_1\\ {\tilde{N}}_2\\ {\tilde{N}}_3\end{array}\right)$

则有：

命题5：1). $D\tilde{N}=e-N$；2). $M_{ij}=\left(\partial \tilde{N}/\partial \xi ,\partial \tilde{N}/\partial \eta \right)diag\left(h_{i,x}^{-1},h_{j,y}^{-1}\right)$；3). $\Vert M_{ij}\Vert =O\left(h^{-1}\right)$；4). $U_{ij}{V}_{ij}={\bar{v}}_{ij}^{\text{T}}N{N}^{\text{T}}{\bar{u}}_{ij}={\bar{v}}_{ij}\left(D\stackrel{⌢}{N}{D}^{\text{T}}+\bar{N}\right){\bar{u}}_{ij}$；5). $\nabla V_{ij}\cdot \nabla U_{ij}={\bar{v}}_{ij}^{\text{T}}D{M}_{ij}{M}_{ij}^{\text{T}}{D}^{\text{T}}{\bar{u}}_{ij}$。

证明：1)，2)，3). 利用形函数的归一性，即 $N_{0,0}+N_{0,1}+N_{1,0}+N_{1,1}=1$[24] ，直接验证即得。

4). 只需验证 $N{N}^{\text{T}}=D\stackrel{⌢}{N}{D}^{\text{T}}+\bar{N}$。

5). ${\nabla N|}_{e_{ij}}=D{M}_{ij}$， $\nabla V_{ij}\cdot \nabla U_{ij}={{\bar{v}}_{ij}^{\text{T}}\nabla N{\left(\nabla N\right)}^{\text{T}}|}_{e_{ij}}{\bar{u}}_{ij}={\bar{v}}_{ij}^{\text{T}}D{M}_{ij}{M}_{ij}^{\text{T}}{D}^{\text{T}}{\bar{u}}_{ij}$。证毕

记

$\begin{array}{l}{Q}_{i,j}={\displaystyle {\iint }_e{q}_{ij}\bar{N}{J}_{ij}d\xi d\eta }+{\displaystyle {\int }_{\partial e}{\sigma }_{ij}{p}_{ij}\bar{N}{L}_{ij}ds}=diag\left(q_{i,j}^{\left(1\right)},q_{i,j}^{\left(2\right)},q_{i,j}^{\left(3\right)},q_{i,j}^{\left(1\right)}\right)\\ P_{i,j}={\displaystyle {\iint }_e\left(p_{ij}{M}_{ij}{M}_{ij}^{\text{T}}+q_{ij}\stackrel{⌢}{N}\right)J_{i,j}d\xi d\eta }+{\displaystyle {\int }_{\partial e}{\sigma }_{ij}{p}_{ij}\stackrel{⌢}{N}{L}_{ij}ds}\\ b_{i,j}={\displaystyle {\iint }_e{f}_{ij}{J}_{ij}Nd\xi d\eta }={\left(b_{i,j}^{\left(1\right)},b_{i,j}^{\left(2\right)},b_{i,j}^{\left(3\right)},b_{i,j}^{\left(4\right)}\right)}^{\text{T}}\end{array}$

则有：

命题6： $\Vert P_{i,j}\Vert =O\left(\Vert J_{ij}{M}_{ij}{M}_{ij}^{\text{T}}\Vert \right)=O\left(1\right)$， $\Vert Q_{i,j}\Vert =O\left(\Vert J_{ij}\Vert \right)=O\left(h^2\right)$。

证明：根据命题4和命题5结论2)，3)可得。证毕

在 $e_{ij}$上，用 $U_{ij},V_{ij}$取代 $u,v$，则得到问题(4)的基于四角剖分的有限元离散格式 [22] ：

${\displaystyle {\iint }_{e_{ij}}\left[p\left(\nabla u\cdot \nabla v\right)+quv\right]dxdy}+{\displaystyle {\int }_{\partial e_{ij}}\sigma puvds}-{\displaystyle {\iint }_{e_{ij}}fvdxdy}={\bar{v}}_{i,j}^{\text{T}}\left[\left(D{P}_{i,j}{D}^{\text{T}}+Q_{i,j}\right){\bar{u}}_{i,j}-b_{i,j}\right]=0$

从而(5)式变成

${\displaystyle \underset{i=1}{\overset{n+1}{\sum }}{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}{\bar{v}}_{i,j}^{\text{T}}}}\left(D{P}_{i,j}{D}^{\text{T}}+Q_{i,j}\right){\bar{u}}_{i,j}={\displaystyle \underset{i=1}{\overset{n+1}{\sum }}{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}{\bar{v}}_{i,j}^{\text{T}}}}{b}_{i,j}$(6)

3.1.2. 基于三角剖分的Lagrange有限元离散格式

将 $e_{ij}$剖分为两个三角形单元：以 $X_{i,j},X_{i,j+1},X_{i+1,j+1}$为顶点的三角形单元记为 $e_{ij}^{\left(1\right)}$，以 $X_{i,j},X_{i+1,j},X_{i+1,j+1}$为顶点的三角形单元记为 $e_{ij}^{\left(2\right)}$，上述变换下 $e^{\left(1\right)}=\left\{\left(\xi ,\eta \right)\in e,\xi \le \eta \right\}$为 $e_{ij}^{\left(1\right)}$的像， $e^{\left(2\right)}=\left\{\left(\xi ,\eta \right)\in e,\xi \ge \eta \right\}$为 $e_{ij}^{\left(2\right)}$的像，在 $e_{ij}^{\left(\alpha \right)}$上取Lagrange形函数 $N_{\beta }^{\left(\alpha \right)}=N_{\beta }^{\left(\alpha \right)}\left({\xi }_1^{\left(\alpha \right)},{\xi }_2^{\left(\alpha \right)},{\xi }_3^{\left(\alpha \right)}\right)$， $\beta =1,2,3$， ${\xi }_1^{\left(\alpha \right)},{\xi }_2^{\left(\alpha \right)},{\xi }_3^{\left(\alpha \right)}$是形函数的自然坐标，则有 ${\xi }_1^{\left(\alpha \right)}+{\xi }_2^{\left(\alpha \right)}+{\xi }_3^{\left(\alpha \right)}=1$， $\alpha =1,2$[23] 。

记 $N^{\left(1\right)}={\left(N_1^{\left(1\right)},0,N_2^{\left(1\right)},N_3^{\left(1\right)}\right)}^{\text{T}}$， $N^{\left(2\right)}={\left(N_1^{\left(2\right)},N_2^{\left(2\right)},0,N_3^{\left(2\right)}\right)}^{\text{T}}$， $\bar{N}=diag\left(N^{\left(1\right)}+N^{\left(2\right)}\right)$， ${\tilde{N}}^{\left(1\right)}={\left(N_1^{\left(1\right)},0,N_2^{\left(1\right)}\right)}^{\text{T}}$， ${\tilde{N}}^{\left(2\right)}={\left(N_1^{\left(2\right)},N_2^{\left(2\right)},0\right)}^{\text{T}}$， $\stackrel{⌢}{N}={\tilde{N}}^{\left(1\right)}{\tilde{N}}^{\left(1\right)\text{T}}+{\tilde{N}}^{\left(2\right)}{\tilde{N}}^{\left(2\right)\text{T}}-diag\left({\tilde{N}}^{\left(1\right)}+{\tilde{N}}^{\left(2\right)}\right)$， ${\bar{u}}_{ij}={\left(u_{i-1,j-1},u_{i-1,j},u_{i,j-1},u_{i,j}\right)}^{\text{T}}$， ${\bar{v}}_{ij}={\left(v_{i-1,j-1},v_{i-1,j},v_{i,j-1},v_{i,j}\right)}^{\text{T}}$，

$U_{ij}={\bar{u}}_{ij}^{\text{T}}\left(N^{\left(1\right)},N^{\left(2\right)}\right)$， $V_{ij}={\bar{v}}_{ij}^{\text{T}}(N^{\left(1\right)},N^{(\; 2\; )})$

$\bar{D}=\frac{1}{2}\left(\begin{array}{ccc}1& 0& 1\\ 1& 1& 0\\ 0& 1& 1\end{array}\right)$， $M_{ij}^{\left(\alpha \right)}=\left(\frac{\partial {\tilde{N}}^{\left(\alpha \right)}}{\partial \xi },\frac{\partial {\tilde{N}}^{\left(\alpha \right)}}{\partial \eta }\right)\left(\begin{array}{cc}{h}_{i,x}^{-1}& 0\\ 0& h_{j,y}^{-1}\end{array}\right)$， ${\bar{M}}_{ij}^{\left(\alpha \right)}={\left(\frac{\partial N^{\left(\alpha \right)}}{\partial x},\frac{\partial N^{\left(\alpha \right)}}{\partial y}\right)|}_{e_{ij}^{\left(\alpha \right)}}$， $\alpha =1,2$

则有：

命题7：1). $D\bar{D}{\tilde{N}}^{\left(\alpha \right)}=e-N^{\left(\alpha \right)}$；2). $\Vert M_{ij}^{\left(\alpha \right)}\Vert =O\left(h^{-1}\right)$；

3). ${\bar{M}}_{ij}^{\left(\alpha \right)}=-D\bar{D}{M}_{ij}^{\left(\alpha \right)}$；4). $U_{ij}{V}_{ij}={\bar{v}}_{ij}^{\text{T}}\left[D\bar{D}\stackrel{⌢}{N}{\bar{D}}^{\text{T}}{D}^{\text{T}}+\bar{N}\right]{\bar{u}}_{ij}$；

5). $\nabla V_{ij}\cdot \nabla U_{ij}={\bar{v}}_{ij}^{\text{T}}D\bar{D}\left(M_{ij}^{\left(1\right)}{M}_{ij}^{\left(1\right)\text{T}}+M_{ij}^{\left(2\right)}{M}_{ij}^{\left(2\right)\text{T}}\right){\bar{D}}^{\text{T}}{D}^{\text{T}}{\bar{u}}_{ij}$；

证明：1)，2)，利用形函数的归一性，即 $N_1^{\left(\alpha \right)}+N_2^{\left(\alpha \right)}+N_3^{\left(\alpha \right)}=1$[24] ，直接验证即得。

3) 根据1)，可验证

${\bar{M}}_{ij}^{\left(\alpha \right)}={\left(\frac{\partial N^{\left(\alpha \right)}}{\partial x},\frac{\partial N^{\left(\alpha \right)}}{\partial y}\right)|}_{e_{ij}^{\left(\alpha \right)}}=-D\bar{D}{\left(\frac{\partial {\tilde{N}}^{\left(\alpha \right)}}{\partial x},\frac{\partial {\tilde{N}}^{\left(\alpha \right)}}{\partial y}\right)|}_{e_{ij}^{\left(\alpha \right)}}=-D\bar{D}{M}_{ij}^{(\; \alpha \; )}$

4). $U_{ij}{V}_{ij}={\bar{v}}_{ij}^{\text{T}}\left(N^{\left(1\right)}{N}^{\left(1\right)\text{T}}+N^{\left(2\right)}{N}^{\left(2\right)\text{T}}\right){\bar{u}}_{ij}$，容易验证 $N^{\left(1\right)}{N}^{\left(1\right)\text{T}}+N^{\left(2\right)}{N}^{\left(2\right)\text{T}}=D\bar{D}\stackrel{⌢}{N}{\bar{D}}^{\text{T}}{D}^{\text{T}}+{\bar{N}}^{\left(1\right)}+{\bar{N}}^{\left(2\right)}$，所以结论成立。

5). $\nabla V_{ij}\cdot \nabla U_{ij}={\bar{v}}_{ij}^{\text{T}}{\left[{\bar{M}}_{ij}^{\left(1\right)}{\bar{M}}_{ij}^{\left(1\right)\text{T}}+{\bar{M}}_{ij}^{\left(2\right)}{\bar{M}}_{ij}^{\left(2\right)\text{T}}\right]|}_{e_{ij}}{\bar{u}}_{ij}$，利用结论3)即可验证。证毕

记

$\begin{array}{l}{\bar{b}}_{i,j}={\displaystyle {\iint }_{e^{\left(\alpha \right)}}{f}_{ij}{J}_{ij}}\left(N^{\left(1\right)}+N^{\left(2\right)}\right)d\xi d\eta ={\left({\bar{b}}_{i,j}^{\left(1\right)},{\bar{b}}_{i,j}^{\left(2\right)},{\bar{b}}_{i,j}^{\left(3\right)},{\bar{b}}_{i,j}^{\left(4\right)}\right)}^{\text{T}}\\ Q_{i,j}^{\left(\alpha \right)}={\displaystyle {\iint }_{e^{\left(\alpha \right)}}{q}_{ij}{\bar{N}}^{\left(\alpha \right)}{J}_{ij}d\xi d\eta }+{\displaystyle {\int }_{\partial e^{\left(\alpha \right)}}{\sigma }_{ij}{p}_{ij}{\bar{N}}^{\left(\alpha \right)}{L}_{ij}ds}\\ P_{i,j}^{\left(\alpha \right)}={\displaystyle {\iint }_{e^{\left(\alpha \right)}}\left(p_{ij}{M}_{ij}^{\left(\alpha \right)}{M}_{ij}^{\left(\alpha \right)\text{T}}+q_{ij}{\stackrel{⌢}{N}}^{\left(\alpha \right)}\right)J_{i,j}d\xi d\eta }+{\displaystyle {\int }_{\partial e^{\left(\alpha \right)}}{\sigma }_{ij}{p}_{ij}{\stackrel{⌢}{N}}^{\left(\alpha \right)}}{L}_{ij}ds,\alpha =1,2\\ {\bar{Q}}_{i,j}=Q_{i,j}^{\left(1\right)}+Q_{i,j}^{\left(2\right)}=diag\left({\bar{q}}_{i,j}^{\left(1\right)},{\bar{q}}_{i,j}^{\left(2\right)},{\bar{q}}_{i,j}^{\left(3\right)},{\bar{q}}_{i,j}^{\left(4\right)}\right)\\ {\bar{P}}_{i,j}=\bar{D}\left[P_{i,j}^{\left(1\right)}+P_{i,j}^{\left(2\right)}\right]{\bar{D}}^{\text{T}}\end{array}$

则有：

命题8： $\Vert {\bar{P}}_{i,j}\Vert =O\left(\Vert J_{ij}\left(M_{ij}^{\left(1\right)}{M}_{ij}^{\left(1\right)\text{T}}+M_{ij}^{\left(2\right)}{M}_{ij}^{\left(2\right)\text{T}}\right)\Vert \right)=O\left(1\right)$， $\Vert {\bar{Q}}_{i,j}\Vert =O\left(\Vert J_{ij}\Vert \right)=O\left(h^2\right)$。

证明：根据命题4和命题7结论2)，3)可得。证毕

在 $e_{ij}$上，用 $U_{ij},V_{ij}$作为 $u,v$的近似，则也可以得到问题(4)的基于三角剖分的有限元离散格式：

${\displaystyle {\iint }_{e_{ij}}\left[p\left(\nabla u\cdot \nabla v\right)+quv\right]dxdy}+{\displaystyle {\int }_{\partial e_{ij}}\sigma puvds}-{\displaystyle {\iint }_{e_{ij}}fvdxdy}={\bar{v}}_{i,j}^{\text{T}}\left[\left(D{\bar{P}}_{i,j}{D}^{\text{T}}+{\bar{Q}}_{i,j}\right){\bar{u}}_{i,j}-{\bar{b}}_{i,j}\right]=0$

从而(5)式变成

${\displaystyle \underset{i=1}{\overset{n+1}{\sum }}{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}{\bar{v}}_{i,j}^{\text{T}}}}\left(D{\bar{P}}_{i,j}{D}^{\text{T}}+{\bar{Q}}_{i,j}\right){\bar{u}}_{i,j}={\displaystyle \underset{i=1}{\overset{n+1}{\sum }}{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}{\bar{v}}_{i,j}^{\text{T}}}}{\bar{b}}_{i,j}$(7)

3.2. 二次Lagrange有限元离散系统

记 $v_i={\left(v_{i,1},\cdots ,v_{i,m}\right)}^{\text{T}}$， ${\tilde{v}}_i={\left(v_{i,0},v_i^{\text{T}},v_{i,m+1}\right)}^{\text{T}}$， $v={\left(v_1^{\text{T}},\cdots ,v_n^{\text{T}}\right)}^{\text{T}}$， $\tilde{v}={\left({\tilde{v}}_0^{\text{T}},\cdots ,{\tilde{v}}_{n+1}^{\text{T}}\right)}^{\text{T}}$

$u_i={\left(u_{i,1},\cdots ,u_{i,m}\right)}^{\text{T}}$， ${\tilde{u}}_i={\left(u_{i,0},u_i^{\text{T}},u_{i,m+1}\right)}^{\text{T}}$， $u={\left(u_1^{\text{T}},\cdots ,u_n^{\text{T}}\right)}^{\text{T}}$， $\tilde{u}={\left({\tilde{u}}_0^{\text{T}},\cdots ,{\tilde{u}}_{n+1}^{\text{T}}\right)}^{\text{T}}$

$Q_i^{\left(s\right)}=diag\left(q_{i,2}^{\left(2s-1\right)}+q_{i,1}^{\left(2s\right)},\cdots ,q_{i,m+1}^{\left(2s-1\right)}+q_{i,m}^{\left(2s\right)}\right)$， ${\tilde{Q}}_i^{\left(s\right)}=diag\left(q_{i,1}^{\left(2s-1\right)},Q_i^{\left(s\right)},q_{i,m+1}^{\left(2s\right)}\right)$， $s=1,2$

$\tilde{Q}=diag\left({\tilde{Q}}_1^{\left(1\right)},{\tilde{Q}}_2^{\left(1\right)}+{\tilde{Q}}_1^{\left(2\right)},\cdots ,{\tilde{Q}}_{n+1}^{\left(1\right)}+{\tilde{Q}}_n^{\left(2\right)},{\tilde{Q}}_{n+1}^{\left(2\right)}\right)$， $Q=diag(Q_2^{\left(1\right)}+Q_1^{\left(2\right)},\cdots ,Q_{n+1}^{\left(1\right)}+Q_n^{(\; 2\; )}$

$b_i^{\left(s\right)}={\left(b_{i,2}^{\left(2s-1\right)}+b_{i,1}^{\left(2s\right)},\cdots ,b_{i,m+1}^{\left(2s-1\right)}+b_{i,m}^{\left(2s\right)}\right)}^{\text{T}}$， ${\tilde{b}}_i^{\left(s\right)}={\left(b_{i,1}^{\left(2s-1\right)},b_i^{\left(s\right)\text{T}},q_{i,m+1}^{\left(2s\right)}\right)}^{\text{T}}$， $s=1,2$

$\tilde{b}={\left({\tilde{b}}_1^{\left(1\right)\text{T}},{\tilde{b}}_2^{\left(1\right)\text{T}}+{\tilde{b}}_1^{\left(2\right)\text{T}},\cdots ,{\tilde{b}}_{n+1}^{\left(1\right)\text{T}}+{\tilde{b}}_n^{\left(2\right)\text{T}},{\bar{b}}_{n+1}^{\left(2\right)\text{T}}\right)}^{\text{T}}$， $b={\left(b_2^{\left(1\right)\text{T}}+b_1^{\left(2\right)\text{T}},\cdots ,b_{n+1}^{\left(1\right)\text{T}}+b_n^{\left(2\right)\text{T}}\right)}^{\text{T}}$

$P_i=diag\left(P_{i,1},\cdots ,P_{i,m+1}\right)$， $P=diag\left(P_1,\cdots ,P_{n+1}\right)$，则有：

命题9：离散格式(6)的矩阵形式为：

${\tilde{v}}^{\text{T}}\left[\left(\bar{Z}P{\bar{Z}}^{\text{T}}+\tilde{Q}\right)\tilde{u}-\tilde{b}\right]=\text{0}$(8)

(8)式中 $\bar{Z}={\left(I_{n+1},0_{n+1}\right)}^{\text{T}}\otimes Z_0+{\left(0_{n+1},I_{n+1}\right)}^{\text{T}}\otimes Z_1$， $Z_i={\left(I_{m+1},0_{m+1}\right)}^{\text{T}}\otimes {\sigma }_{2i}+{\left(0_{m+1},I_{m+1}\right)}^{\text{T}}\otimes {\sigma }_{2i+1}$， $i=0,1$。

证明：容易验证：

${\displaystyle \underset{i=1}{\overset{n+1}{\sum }}{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}{\bar{v}}_{ij}^{\text{T}}D{P}_{ij}{D}^{\text{T}}{\bar{u}}_{ij}}}={\displaystyle \underset{i=1}{\overset{n+1}{\sum }}\left({\tilde{v}}_{i-1}^{\text{T}},{\tilde{v}}_i^{\text{T}}\right)\left(\begin{array}{c}{\bar{Z}}_0\\ {\bar{Z}}_1\end{array}\right)P_i\left({\bar{Z}}_0^{\text{T}},{\bar{Z}}_1^{\text{T}}\right)\left(\begin{array}{c}{\tilde{u}}_{i-1}\\ {\tilde{u}}_i\end{array}\right)}={\tilde{v}}^{\text{T}}\bar{Z}P{\bar{Z}}^{\text{T}}\tilde{u}$

${\displaystyle \underset{i=1}{\overset{n+1}{\sum }}{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}{\bar{v}}_{ij}^{\text{T}}{Q}_{i,j}{\bar{u}}_{ij}}}={\displaystyle \underset{i=1}{\overset{n+1}{\sum }}\left({\tilde{v}}_{i-1}^{\text{T}}{\tilde{Q}}_i^{\left(1\right)}{\tilde{u}}_{i-1}+{\tilde{v}}_i^{\text{T}}{\tilde{Q}}_i^{\left(2\right)}{\tilde{u}}_i\right)}={\tilde{v}}^{\text{T}}\tilde{Q}\tilde{u}$

${\displaystyle \underset{i=1}{\overset{n+1}{\sum }}{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}{\bar{v}}_{ij}^{\text{T}}{b}_{ij}}}={\displaystyle \underset{i=1}{\overset{n+1}{\sum }}\left({\tilde{v}}_{i-1}^{\text{T}}{\tilde{b}}_i^{\left(1\right)}+{\tilde{v}}_i^{\text{T}}{\tilde{b}}_i^{\left(2\right)}\right)}={\tilde{v}}^{\text{T}}\tilde{b}$

从而将问题(4)的离散格式(6)写成公式(8)。证毕

根据式(8)，有

$\left(\bar{Z}P{\bar{Z}}^{\text{T}}+\tilde{Q}\right)\tilde{u}=\tilde{b}$(9)

将(9)式中与边界网格点的值 $u_{0,j},u_{n+1,j},u_{i,0},u_{i,m+1}$相关的部分移到右端，就得到问题(4)基于四角剖分的有限元离散系统：

$Au=b$(10)

式(10)中，刚度矩阵为

$A=ZP{Z}^{\text{T}}+Q$(11)

类似地，可以根据离散格式(7)得到问题(4)基于三角剖分的有限元离散系统：

$\bar{A}u=\left(Z\bar{P}{Z}^{\text{T}}+\bar{Q}\right)u=\bar{b}$(12)

式(12)中 ${\bar{P}}_i=diag\left({\bar{P}}_{i,1},\cdots ,{\bar{P}}_{i,m+1}\right)$， $\bar{P}=diag\left({\bar{P}}_1,\cdots ,{\bar{P}}_{n+1}\right)$

${\bar{Q}}_i^{\left(s\right)}=diag\left({\bar{q}}_{i,2}^{\left(2s-1\right)}+{\bar{q}}_{i,1}^{\left(2s\right)},\cdots ,{\bar{q}}_{i,m+1}^{\left(2s-1\right)}+{\bar{q}}_{i,m}^{\left(2s\right)}\right)$， ${\bar{b}}_i^{\left(s\right)}={\left({\bar{b}}_{i,2}^{\left(2s-1\right)}+{\bar{b}}_{i,1}^{\left(2s\right)},\cdots ,{\bar{b}}_{i,m+1}^{\left(2s-1\right)}+{\bar{b}}_{i,m}^{\left(2s\right)}\right)}^{\text{T}}$， $s=1,2$

$\bar{Q}=diag\left({\bar{Q}}_2^{\left(1\right)}+{\bar{Q}}_1^{\left(2\right)},\cdots ,{\bar{Q}}_{n+1}^{\left(1\right)}+{\bar{Q}}_n^{\left(2\right)}\right)$， $\bar{b}={\left({\bar{b}}_2^{\left(1\right)\text{T}}+{\bar{b}}_1^{\left(2\right)\text{T}},\cdots ,{\bar{b}}_{n+1}^{\left(1\right)\text{T}}+{\bar{b}}_n^{\left(2\right)\text{T}}\right)}^{\text{T}}$

不论是三角剖分，还是四角剖分，离散系统的总刚度矩阵都具有式(1)的形式，这个形式有利于评估条件数，便于设计预条件子和算法。

4. 对有限元离散系统的病因抑制方法及其预处理效果

根据命题6、8，(10)和(12)式中，分别有 $\Vert ZP{Z}^{\text{T}}\Vert \gg \Vert Q\Vert$和 $\Vert Z\bar{P}{Z}^{\text{T}}\Vert \gg \Vert \bar{Q}\Vert$，所以在总刚度矩阵 $A$和 $\bar{A}$中， $ZP{Z}^{\text{T}}$和 $Z\bar{P}{Z}^{\text{T}}$分别是 $A$与 $\bar{A}$的主体， $Cond\left(A\right)\approx Cond\left(ZP{Z}^{\text{T}}\right)$， $Cond\left(\bar{A}\right)\approx Cond\left(Z\bar{P}{Z}^{\text{T}}\right)$；矩阵 $P$和 $\bar{P}$的元素由 $p=p\left(x,y\right)$的振幅，以及网格大小决定，几乎不受 $m,n$影响，所以 $Cond\left(P\right)$和 $Cond\left(\bar{P}\right)$关于 $A$与 $\bar{A}$的阶数都是一致有界的，根据定义1、3和命题2，(10)和(12)式分别说明了刚度矩阵 $A$与 $\bar{A}$都具有(1)式的病态结构。

(10)和(12)式中， $Q$和 $\bar{Q}$对刚度矩阵的病态影响很小。根据命题3，由病态因子 $Z$表达的原发病态是离散系统的主要病态， $P,Q$和 $\bar{P},\bar{Q}$等表达的继发病态是散系统的次要病态，原发病态是微分算子造成的，是本质的，离散精度越高，刚度矩阵的阶数越大，刚度矩阵条件数就越大，原发病态越严重；继发病态来自数值方法和问题(4)本身，是非本质的，可以随着应用问题的不同而不同，可以在应用中调整。

综上所述，(10)和(12)说明，基于不同剖分的离散系统，有类似的病态结构和相同的病态因子 $Z$，即有相同的原发病态。

使用 $Z$的去病因子 $H$作为预条件子，对方程(10)、(12)进行预处理，其中方程(10)化成

$H^{-1}A{H}^{-\text{T}}\tilde{u}=\left(H^{-1}ZP{Z}^{\text{T}}{H}^{-\text{T}}+H^{-1}Q{H}^{-\text{T}}\right)\tilde{u}=\tilde{b}$(13)

其中 $\tilde{u}=H^{\text{T}}u$， $\tilde{b}=H^{-1}b$。

根据命题1的结论2)，有 $Cond\left(H^{-1}A{H}^{-\text{T}}\right)\approx Cond\left(H^{-1}ZP{Z}^{\text{T}}{H}^{-\text{T}}\right)\le Cond\left(P\right)$，因此， $H{H}^{\text{T}}$是最优预条件子 [5] 。

使用 $Z$的去病因子 $H$作为预条件子，对方程(12)进行预处理，也有类似的结果。

$H,H^{-1}$有简单确定的结构，都是对角矩阵与离散正弦变换矩阵的直积的积，因此预条件子的构造以及预处理计算不会明显增加计算成本，其中预处理需要的计算是矩阵 $H^{\text{T}},H^{-1}$与向量的乘积，相当于对向量实施快速离散正弦变换，每次需要的计算操作数为 $O\left(nm{\log }_2\left(nm\right)\right)$，所以每个迭代步的计算量没有显著增加。预处理后，方程(13)的系数矩阵仍然是正定对称矩阵，所以仍然可以与经典Krylov子空间方法等求解方法相结合，形成预处理求解方法 [4] 。对方程(12)进行同样的预处理，也有类似的结果。

5. 例证

对于积分形式的二维泊松方程边值问题(4) [22] ，如果 $\sigma =q=0$，p为常数，取2次Lagrange形函数为： $N_1^{\left(1\right)}=1-\eta$， $N_2^{\left(1\right)}=\xi$， $N_3^{\left(1\right)}=\eta -\xi$， $N_1^{\left(2\right)}=1-\xi$， $N_2^{\left(2\right)}=\xi -\eta$， $N_3^{\left(2\right)}=\eta$，并取 $h_x=h_{i,x}$， $h_y=h_{j,y}$，记 $\bar{h}=h_y{h}_x^{-1}$， $B_n={\left(-1,2,-1\right)}_n$， $C_n={\left(0,1,1\right)}_n$，则

$P_{i,j}=\frac{p}{8}\left(\begin{array}{ccc}2\bar{h}+{\bar{h}}^{-1}& {\bar{h}}^{-1}& \\ {\bar{h}}^{-1}& 2{\bar{h}}^{-1}& {\bar{h}}^{-1}\\ & {\bar{h}}^{-1}& 2\bar{h}+{\bar{h}}^{-1}\end{array}\right)$

$P=\frac{p}{8}{I}_{mn}\otimes \left(\begin{array}{ccc}2\bar{h}+{\bar{h}}^{-1}& {\bar{h}}^{-1}& \\ {\bar{h}}^{-1}& 2{\bar{h}}^{-1}& {\bar{h}}^{-1}\\ & {\bar{h}}^{-1}& 2\bar{h}+{\bar{h}}^{-1}\end{array}\right)$

$A=ZP{Z}^{\text{T}}=\frac{p}{2}\left[2\bar{h}\left(\begin{array}{ccc}{B}_n& & \\ & \ddots & \\ & & B_n\end{array}\right)+{\bar{h}}^{-1}\left(\begin{array}{cccc}4I_n& -C_n& & \\ -C_n^{\text{T}}& \ddots & \ddots & \\ & \ddots & \ddots & -C_n\\ & & -C_n^{\text{T}}& 4I_n\end{array}\right)\right]$

$Cond\left(ZP{Z}^{\text{T}}\right)\le Cond\left(P\right)Cond\left(Z{Z}^{\text{T}}\right)$

$Cond\left(H^{-1}A{H}^{-\text{T}}\right)\le Cond\left(P\right)\le 4\max \left(\bar{h},{\bar{h}}^{-1}\right)/\min \left(\bar{h},{\bar{h}}^{-1}\right)$，即预处理后，矩阵的条件数一致有界。

特别地，如果取 $h_x=h_y$，则 $\bar{h}=1$， $Cond\left(H^{-1}A{H}^{-\text{T}}\right)\le 4$，即预处理后，矩阵的条件数小于4，因此，去病因子 $H$是最优的预条件子。

6. 结论

1) 二维泊松方程边值问题基于三角、四角剖分的Lagrange有限元离散系统虽然不同，其系数矩阵却有类似的病态结构，有相同的病态因子，它的病态由原发病态和继发病态组成，原发性病态是由微分算子引起的，由病态因子表达，是主要的、本质的，不受剖分方法影响，可以使用相同的去病因子进行预处理；继发病态是由数值方法和边值问题的系数函数引起，是次要的、非本质的，可以在应用中调整。

2) 去病因子是最优预条件子 [5] ，精准针对病态因子起作用，可消除原发性病态，几乎将矩阵的条件数降为常数，相当于精准地抑制了病态因子的致病作用。去病因子的特别结构，使得预处理过程基本不增加计算操作数，并且预处理后矩阵保持正定对称，有利于与其他求解方法结合。

3) 病态因子及其去病因子是相对独立的，与边值问题中函数p，q，f等没有关系，不论三角剖分，还是四角剖分，都不影响病态因子和去病因子，但是剖分方法、步长可以影响继发病态。

4) 去病因子的设计是针对病态因子(无须针对整个刚度矩阵)，只要从刚度矩阵中分离出病态因子，就可以得到去病因子，因此去病因子的获取难度有限。

综上所述，与目前主流的预处理方法相比，本文提出的病因抑制方法，具有难度低，不增加计算成本，精准抑制原发病态等特点。

基金项目

国家自然基金资助项目(61473329, 61601125)。