特征值问题虚拟元方法发展综述

孟 健, 梅立泉

(西安交通大学数学与统计学院，西安 710049)

0 引言

偏微分方程特征值问题是非常重要的研究课题，在科学和工程应用领域有重要意义，例如，非线性方程谱理论、结构力学、流体稳定性、电子结构计算、控制理论和逆散射问题[1–7]。特征值问题的数值方法和数学理论涉及十分广泛，迄今为止已有一些有重要意义的著作，如Chatelin 研究了Banach 空间中线性算子特征值的数值逼近，把泛函分析谱理论框架应用到特征值问题有限元方法[8]；Babˇuska 和Osborn 进一步对特征值问题协调有限元方法和混合有限元方法进行了系统的总结[9]；文献[10]讨论了特征值问题有限元方法的渐进展开和外推；最新的一些著作，可参考文献[11—13]。而且已有许多数值方法用来求解特征值问题，如有限差分方法[14]、有限元方法[9,11–12]、谱方法[15]、杂交高阶方法[16]、间断Galerkin 方法[17–18]和弱Galerkin 方法[19]等。

文献[20]结合有限元和模拟有限差分方法[21–24]引入了虚拟元方法，阐述了虚拟元方法的主要思想：

1) 局部离散空间包含多项式函数和由局部微分问题决定的非多项式函数(“虚拟函数”)；

2) 为了得到在一般多边形网格上完全可计算的虚拟元格式，虚拟元方法需小心选取自由度。

作为一种新型偏微分方程数值方法，虚拟元方法具有如下特点：

1) 允许一般多边形网格，有好的网格适应性；

2) 易构造高阶元离散空间。

相较于传统有限元方法，虚拟元法有好的网格适应性和数值稳定性，降低了复杂几何区域上网格生成和自适应算法的难度。如果采用三角形网格，虚拟元方法也等价于传统有限元方法。除此之外，虚拟元格式与有限元格式有一定的关联性，所以虚拟元方法理论分析可以借鉴有限元方法。除此之外，高阶/维虚拟元空间的构造与低阶/维的构造方式统一，其程序编写也相对容易。近年来，各类虚拟元空间已经被发展来处理各种数学与工程问题，例如，非协调虚拟元[25–26]、H(div)和H(curl)虚拟元[27–29]、Cm虚拟元[30–34]。并且，许多学者研究了特征值问题的虚拟元方法，包括Steklov 特征值问题[35–37]、Laplacian 特征值问题[38–42]、声音振动问题[43]、Kirchhoff板振动和屈曲问题[33,44–45]、透射特征值问题[32]和弹性特征值问题[46]。

本文将总结特征值问题虚拟元方法的主要文献，简要说明虚拟元方法求解特征值问题的特点，提出今后的研究方向和建议，促进特征值问题数值方法的研究。

本文结构如下：第1 节引入虚拟元方法，定义所需虚拟元空间和自由度，包括协调虚拟元空间、H(div)虚拟元空间和C1虚拟元空间；第2 节介绍特征值问题虚拟元方法的研究现状；第3 节是本文的总结与展望。

1 虚拟元方法介绍

有限元方法的出现对当代计算数学有重大意义，在水坝和飞机的设计、天体物理和流体力学等问题有重要应用，有限元方法常用的网格单元有三角形和四边形。近年来，有限元方法的推广也是数学和工程计算中的热点问题，为了处理复杂区域，许多学者把有限元方法推广到一般的多边形网格情形[20,47–49]。多边形网格容易做网格自适应，可以用来提升数值方法的精度和处理结构复杂的解，例如，有奇异性的解和界面问题，设计多边形单元数值方法面临的关键问题是如何构造离散函数空间。在传统协调有限元方法中，离散函数空间包含在偏微分方程真解函数空间中，当使用多边形网格时，会导致其上的基函数求解困难，为了满足空间的协调性，基函数的形式并不一定是多项式[50]。

结合模拟有限差分方法和有限元方法，虚拟元方法被引入[20]，如上讨论，虚拟元空间包含多项式函数和局部适定微分问题决定的非多项式函数。在虚拟元方法中，实际不需要计算这些非多项式函数，所以称作“虚拟函数”，可以通过考虑虚拟元空间的定义和自由度来使用这些“虚拟函数”。对于构造虚拟元离散双线性形式，可以通过一致项和稳定项来保证最后得到矩阵问题的适定性。根据作者理解，虚拟元方法的优势主要有：

1) 相较于传统有限元方法，降低了复杂几何区域上网格生成和自适应算法难度；

2) 高维和高正则性Cm虚拟元空间可以由低维和低正则性虚拟元空间递归定义；

3) 容易构造满足物理性质的空间和数值格式，如满足无散度性质[51–52]等。

下面引入特征值问题使用的三类虚拟元空间。

1.1 协调虚拟元空间

令Ω ⊆R2是一个有界区域，边界∂Ω上的单位外法向记为ν。令Th是Ω上多边形网格剖分，设E是其中某个单元，Eh是剖分单元所有边e的集合。记he、hE、|E|分别为边e的长度、E的直径和E的面积，且网格尺寸h= max{hE:E ∈Th}。假设网格Th满足如下正则性条件，存在一致正常数σ，使得对每个单元E，有：

(A1) 单元E是星形的；

(A2) 对于每条边e ∈∂E，有he ≥σhE。

由文献[20,53]，局部初始空间定义如下

对任意v ∈～VE，定义如下线性算子集合：

LB1:LB1(v)表示v在E所有顶点处的函数值；

LB2:LB2(v)表示v在E的每条边内部k −1 个Gauss-Lobatto 积分点处的函数值；LB3:LB3(v)表示v直到k −2 阶矩的值。

1.2 H(div)虚拟元空间

设给定区域Ω的网格剖分Th，满足1.1 节的正则性假设(A1)和(A2)。在文献[28]中，H(div)虚拟元空间的定义被引入，一些基本的推广见文献[27,29,55]。H(div)虚拟元被应用于各种数学问题，例如，Stokes 问题[56]、拟牛顿Stokes 问题[57]、线性和非线性Brinkman 问题[58–59]。H(div)虚拟元空间定义如下

自由度定义为：

(C1) ∫

e σh·npds, ∀p ∈Mk(e), ∀e ∈Eh；

(C2) ∫

E σh·∇pdx, ∀p ∈Mk−1(E){1}, ∀E ∈Th；

(C3) ∫

Erotσhpdx, ∀p ∈Mk−1(E), ∀E ∈Th；

其中Mk(e)表示边e(设中点为xe)上的单位化单项式集合Mk−1(E)表示单元E(设重心为xE)上的单位化单项式集合

1.3 C1 虚拟元空间

众所周知，全局C1协调有限元的构造是非常困难的[60]。由于虚拟元方法的灵活性，C1虚拟元构造变得相对容易。下面阐述C1虚拟元空间的构造方法，这个虚拟元空间已经应用到很多问题[32–34,45]。局部C1初始空间定义如下

LD1:LD1表示v在E所有顶点处的函数值；

LD2:LD2表示∇v在E所有顶点处的函数值。

2 特征值问题虚拟元方法研究现状

本节介绍特征值问题虚拟元方法(协调虚拟元方法、H(div)虚拟元方法、C1虚拟元方法)已有的研究工作，以下按不同的虚拟元方法进行分类。

2.1 协调虚拟元方法

文献[36]讨论了对称Steklov 特征值问题的虚拟元方法，这也是首篇使用虚拟元方法研究特征值问题的文献。在网格正则性假设下，给出了虚拟元格式的谱逼近性质，证明了特征函数和特征值的最优误差估计。因为Steklov 问题是边界特征值问题，也考虑了在边界上特征函数的高阶收敛性。2017 年，Mora 等人[37]利用虚拟元方法网格灵活性考虑了对称Steklov 特征值问题协调虚拟元格式的后验误差估计。定义了一个可计算的虚拟元后验误差估计子，并把它应用于自适应网格剖分。数值实验验证了对称Steklov 特征值问题协调虚拟元格式数值解的最优收敛性和后验误差估计子的有效性，这是第一篇特征值问题虚拟元方法的后验误差分析。之后，Gardini 和Vacca[41]讨论了Laplacian 特征值问题的协调虚拟元方法，给出了两种协调虚拟元格式，区别在于右端项有无稳定项。应用紧算子和非紧算子谱理论分别证明了两种协调虚拟元格式的谱逼近性质和误差估计，在数值实验中，作者说明对于高阶协调虚拟元方法，应用对角稳定项的数值格式收敛性更好，为以后讨论更复杂的特征值问题虚拟元方法奠定了研究基础。应用这篇文章的理论框架，文献[40]考虑了Laplacian 特征问题的非协调虚拟元方法，文献[38]分析了椭圆特征值问题的p与hp虚拟元方法，ˇCert´ık 等人[39]研究了有势函数项的椭圆特征值问题，在数值实验中模拟了与量子振荡相关的Schr¨odinger 问题，验证了虚拟元方法的有效性。最近，Boffi等人[61]讨论了虚拟元方法中稳定化参量对Laplacian 特征值的影响。Meng 和Mei[44]应用协调虚拟元方法讨论了简支Kirchhoff板屈曲特征值问题，给出了虚拟元格式并使用“第一Strang 引理”证明了特征模型的误差分析和最优收敛阶，数值实验验证了理论结果。Mora 和Rivera[46]讨论了弹性特征值问题虚拟元方法的先验和后验误差估计。

2.2 H(div)虚拟元方法

迄今为止，特征值问题H(div)虚拟元方法的参考文献十分有限，第一篇参考文献是使用H(div)虚拟元处理声音振动问题[43]。在网格正则性假设下，证明了虚拟元格式的正确谱逼近和最优误差估计，数值实验验证了理论结果。最近，Meng 等人[42]应用混合虚拟元方法讨论了Laplacian 特征值问题。在Brezzi-Babuska 理论框架下，验证了解空间的弱逼近性、强逼近性和Fortin 条件，从而证明了数值格式的收敛性和最优性，数值实验验证了混合虚拟元方法的最优收敛性。在这篇文章的基础上，Lepe 和Rivera[62]证明了Laplacian 特征值问题混合虚拟元方法的误差估计。

2.3 C1 虚拟元方法

Mora 等人[33]考虑了夹紧Kirchhoff板振动特征值问题，应用非紧算子的谱理论证明了谱逼近和误差估计，数值实验验证了在凸和非凸区域上特征值问题虚拟元方法的收敛性。之后，Mora 和Vel´asquez[32]应用C1虚拟元方法分析了非对称透射特征值问题，这个特征值问题是四阶且非对称的，其理论和数值研究是具有挑战性的研究课题，诸多学者对透射特征值问题的数值方法进行了研究。特别地，C1Argyris 有限元方法已经被用来求解这个特征值问题[63–64]，但全局C1协调有限元的构造是非常困难的[60]。由于虚拟元方法的灵活性，C1虚拟元构造变得相对容易。最后，Mora 和Vel´asquez[45]推广了文献[32—33]中C1虚拟元空间到任意阶(k ≥2)的情况，并处理了夹紧Kirchhoff板屈曲特征值问题。

3 总结与展望

自从虚拟元方法诞生以来，特征值问题虚拟元方法的理论和数值方法取得了一些进展，使用各种虚拟元方法处理了各类特征值问题，但还有很多关键问题亟需解决。

1) 非对称特征值问题的虚拟元方法。非对称特征值问题有重要的应用背景，例如，环境问题和流体力学里中的对流扩散，但非对称特征值问题必须考虑复特征值且无法保证特征值的陡度为1，从现有文献看，其数值方法研究工作略滞后于对称特征值问题。由于工程问题的实际需求，非对称特征值问题虚拟元方法的理论和数值分析有重要意义。例如，来源于非均匀介质逆散射理论的非对称Steklov 特征值问题，特征值能通过远场数据得到并且可以用来分析散射物体的属性，在很多领域有着广泛应用。目前，非对称Steklov 特征值问题数值方法的参考文献十分有限，如有限元方法、非协调Crouzeix-Raviart 元方法、多网格方法和间断Galerkin 方法。其虚拟元方法可以使用已有文献的虚拟元离散空间，但非对称性阻碍了数值方法的理论分析和对复数特征值的求解，这也是非对称特征值问题面临的困难。

2) 高效准确的数值算法。目前，已有部分特征值问题虚拟元方法的文献，下一目的是提出能提高计算效率和减少计算消耗的虚拟元格式。首先，在特征值问题限元方法中，存在很多加速算法，例如，多重网格算法、两空间加速算法和位移反幂算法，都可以推广到虚拟元方法中，并给出相应的理论和数值分析。其次，对于高阶虚拟元方法的构造，应注意虚拟元格式的定义，原始虚拟元格式在使用高阶(k ≥3)虚拟元空间时，数值结果的收敛阶是次优的。所以为构造高效准确的数值方法，也应特别注意虚拟元格式的构造。最后，也可以从问题本身出发，考虑不同的虚拟元方法。

以方形区域Ω= (0,1)2为例，考虑如图1 的三角形网格Th、方形网格Sh、泰森多边形网格Vh和无结构六边形网格Nh上两种方法的自由度个数，表1 给出了在四种网格上两种不同方法的自由度个数，发现以往文献中的方法使用了更少的自由度，提高了计算效率。但随之而来的问题是如何进行理论分析？在以往文献中，作者使用“第一Strang 引理”进行理论分析，但这时不能直接使用这个引理，需要考虑另外的理论分析方法。最近的文献给出了这个问题最低阶虚拟元方法谱逼近理论，并在如图2 三维网格上考虑了三维透射特征值，但因为问题的非自共轭性，文中还未证明虚拟元方法的误差估计。

表1 不同方法的自由度个数

图1 二维网格描述

图2 三维网格描述

3) 谱理论的推广。现有特征值问题虚拟元方法所应用的谱理论是紧算子和非紧算子的谱理论框架，也可以考虑把其他算子的谱理论推广到虚拟元方法来处理更复杂的特征值问题，可以推广Fredholm 算子谱理论来处理各向异性的透射特征值问题。

4) 其他特征值问题的虚拟元方法。总结来看，使用虚拟元方法处理的特征值问题包括：Steklov 特征值问题、Laplacian 特征值问题、声音振动问题、Kirchhoff板振动和屈曲问题、透射特征值问题和弹性特征值问题。目前用H(div)虚拟元方法处理特征值问题的文献仅有三篇，但H(div)元在特征值问题中有着重要应用，对于流体结构振动问题，基于对偶形式的标准有限元方法导出的广义特征值问题会引入错误的特征模型，H(div)有限元方法就可以避免这种谱污染，现有的H(div)虚拟元研究为以后讨论更复杂流体结构特征值问题奠定了研究基础。另外，对于三维特征值问题的虚拟元方法数值结果较少，限制了在实际中的应用。所以可将虚拟元方法应用于其他更具实际意义的特征值问题，如流体、输运、电磁场等特征值问题，解决更多实际问题。