求解流固振动Laplace模型的基于Rayleigh商移位反迭代的多网格方案

Multiscale Discretization Scheme Based on the Rayleigh Quotient Inverse Iteration Method for the Laplace Model for Fluid-Solid Vibrations

摘要: 本文讨论流固振动Laplace模型，首先建立该问题的基于Rayleigh商移位反迭代的多网格离散方案，利用该方案将在细网格上求解特征值问题归结为在粗网格上解特征值问题和在细网格上解一系列线性代数系统。其次分析多网格离散方案的误差。最后给出数值算例验证离散方案的有效性。

Abstract: In this paper, for the Laplace model for fluid-solid vibrations we establish the multiscale discretization scheme based on the Rayleigh quotient inverse iteration method. With this scheme, the solution of an eigenvalue problem on fine meshes is reduced to the solution of the eigenvalue problem on coarse meshes and the solutions of linear algebraic systems. Then we analyze and give the error estimates of the proposed scheme. Finally we present numerical experiments to validate the efficiency of the scheme.

1. 引言

1992年许进超首次引入二网格离散技巧并用于求解非对称和非椭圆线性问题(参见文献 [1] [2] )，2001年许进超和周爱辉又将这一技巧应用到了特征值问题(参见文献 [3] )。由于二网格离散技巧能够显著提高近似解精度和计算效率，节省计算机存贮，因此二网格离散及在其基础上发展起来的多网格离散被广泛用于求解特征值问题，如重调和特征值问题(参见 [4] )，半线性特征值问题(参见 [5] )，量子特征值问题(参见 [6] )，Stokes特征值问题(参见 [7] [8] )，Maxwell特征值问题(参见 [9] )，2m阶椭圆特征值问题(参见 [10] )等。在文献 [1] [2] 建立的二网格离散基础上还发展了基于移位反迭代的二网格、多网格离散以及多水平校正方案。基于移位反迭代的二网格离散、多网格离散可看作解矩阵特征值问题的移位反迭代法与有限元方法的结合(参见 [11] [12] [13] )。之后，文献 [14] 对Laplace特征值问题研究了协调有限元基于移位反迭代的多网格方法； [15] [16] 将该方法用于Maxwell特征值问题的棱元逼近； [17] 将这一工作应用到Stokes特征值问题； [18] 对Steklov特征值问题建立了协调有限元多网格方案； [19] 进一步发展该工作，对一般的自共轭特征值问题建立了基于多网格离散的移位反迭代，并用于积分算子特征值问题； [20] 对重调和特征值问题建立了基于移位反迭代的二网格、多网格方案； [21] 对Laplace特征值问题提出了多水平校正方案。

流固振动Laplace模型问题源于一束浸没在不可压缩流体中的平行管的振动模型的计算，在核工程中有着重要作用(参见文献 [22] [23] [24] )。近年来，该问题的数值方法引起了学者们的研究兴趣，如Armentano等 [25] 研究了该问题的hp协调有限元法及自适应算法，给出了先验误差估计和后验误差估计；张宇等 [26] 研究了非协调有限元特征值的可保证下界。在上述工作基础上，本文研究流固振动Laplace模型的基于Rayleigh商移位反迭代的多网格离散方案。利用该方案可以将在细网格上特征值问题的求解归结为在粗网格上解特征值问题和在细网格上解一系列线性代数组。理论分析表明，当粗网格直径H和细

网格直径h满足 $H=O\left(\sqrt[3]{h}\right)$ 时，多网格方案求得的近似特征值能达到渐进最优收敛阶。

本文剩余部分安排如下。在第2节，我们给出本文所需的预备知识。在第3节，我们建立基于Rayleigh商移位反迭代的多网格方案，并对方案进行数值分析。在第4节给出一些数值实验验证方案的有效性。

文中我们用字母C表示与网格尺寸无关的正常数，它在不同的地方表示的值可能不同。为方便起见，用 $a\lesssim b$ 表示 $a\le Cb$。

2. 预备知识

在合理假设下，流固振动Laplace模型可由下述系统描述，其中每个管被模拟为谐波振荡器且硬度为 $\nu$，质量为m，流体视为完全不可压缩的且密度为 $\rho$ ：求 $\omega >0$ (振动频率)， $u\ne 0$ (流体压力)使得

$\begin{array}{r}\hfill \{\begin{array}{l}\Delta u=0在\Omega 中,\hfill \\ \frac{\partial u}{\partial n}=0在{\Gamma }_0上,\hfill \\ \frac{\partial u}{\partial n}=\frac{\rho {\omega }^2}{\nu -m{\omega }^2}\left({\displaystyle {\int }_{{\Gamma }_i}un\text{d}s}\right)\cdot n在{\Gamma }_i上,i=1,2,\cdots ,K,\hfill \end{array}\end{array}$ (2.1)

其中 $\Omega \subset {ℝ}^2$ 是由流体占据的有界多边形区域， ${\Gamma }_0$ 表示 $\Omega$ 的外边界， ${\Gamma }_i\left(i=1,\cdots ,K\right)$ 表示每个管与流体间的交界面， $n={\left(n_1,n_2\right)}^{\text{T}}$ 是 $\Omega$ 边界上的单位外法向量(参见图1)。

Figure 1. Sketch of the two-dimensional domain

图1. 二维区域模型草图

令 $H^t\left(\Omega \right)$ 和 $H^t\left(\partial \Omega \right)$ 分别表示 $\Omega$ 和 $\partial \Omega$ 上t阶Sobolev空间，其上半范数分别记为 ${\left|\cdot \right|}_{t,\Omega }$ 和 ${\left|\cdot \right|}_{t,\partial \Omega }$， $H^0\left(\partial \Omega \right)=L_2\left(\partial \Omega \right)$。令 $\lambda :=\rho {\omega }^2/\left(\nu -m{\omega }^2\right)$， $V:=H^1\left(\Omega \right)/ℝ$。(2.1)的弱形式为：求 $\lambda \in ℝ$， $u\in V$ 满足

$\{\begin{array}{l}a\left(u,v\right)=\lambda b\left(u,v\right),\forall v\in V,\hfill \\ b\left(u,u\right)=1,\hfill \end{array}$ (2.2)

其中 $a\left(u,v\right):={\displaystyle {\int }_{\Omega }\nabla u\cdot \nabla v}$， $b\left(u,v\right):={\displaystyle {\sum }_{i=1}^K\left({\displaystyle {\int }_{{\Gamma }_i}un}\right)\cdot \left({\displaystyle {\int }_{{\Gamma }_i}vn}\right)}$。

显然 $a\left(\cdot ,\cdot \right)$ 是V上连续、对称、椭圆的双线性形式， $b\left(\cdot ,\cdot \right)$ 是非负的。在V上 $H^1$ 半范数 ${\left|\cdot \right|}_{1,\Omega }=\sqrt{a\left(\cdot ,\cdot \right)}$ 等价于 $H^1$ 范数 ${\Vert \cdot \Vert }_{1,\Omega }$， $\sqrt{b\left(\cdot ,\cdot \right)}$ 可以作为V上的半范数。

由文献 [25] [27] 可知，(2.2)的解是由2K个特征对 $\left({\lambda }_j,u_j\right)$ 给出的序列，特征值 ${\lambda }_j$ 均为正，假定特征值按递增排序： $0<{\lambda }_1\le \cdots \le {\lambda }_{2K}$。每个特征值都对应一个特征函数 $u_j\in V$，使得 $\left\{u_1,\cdots ,u_{2K}\right\}$ 是一个线性无关的集合。

令 ${\pi }_h$ 是 $\Omega$ 的一族正规三角形剖分， $h_T$ 为三角形 $T\in {\pi }_h$ 的直径，网格直径 $h=\underset{T\in {\pi }_h}{\max }{h}_T$。令 ${\epsilon }_h\left({\Gamma }_i\right)=\left\{e\right\}$ 表示 ${\pi }_h$ 中位于边界 ${\Gamma }_i$ 上的所有边构成的集合， ${\epsilon }_h\left(\Gamma \right)={\epsilon }_h\left({\Gamma }_1\right)\cup {\epsilon }_h\left({\Gamma }_2\right)\cup \cdots \cup {\epsilon }_h\left({\Gamma }_K\right)$。令 $S_h\subset H^1\left(\Omega \right)$ 是定义在 ${\pi }_h$ 上的分片p次( $p\ge 1$ )多项式空间， $V_h=S_h/ℝ$。注意到对任一 $v\in V+V_h$，由柯西–施瓦茨不等式、迹定理和商空间等价模定理，有

$\begin{array}{c}b\left(v,v\right)={\displaystyle \underset{i=1}{\overset{K}{\sum }}\left({\displaystyle \underset{e\in {\epsilon }_h\left({\text{Γ}}_i\right)}{\sum }{\displaystyle \underset{e}{\int }vn}}\right)\left({\displaystyle \underset{e\in {\epsilon }_h\left({\text{Γ}}_i\right)}{\sum }{\displaystyle \underset{e}{\int }vn}}\right)}={\displaystyle \underset{i=1}{\overset{K}{\sum }}\left\{{\left({\displaystyle \underset{e\in {\epsilon }_h\left({\text{Γ}}_i\right)}{\sum }{\displaystyle \underset{e}{\int }v{n}_1}}\right)}^2+{\left({\displaystyle \underset{e\in {\epsilon }_h\left({\text{Γ}}_i\right)}{\sum }{\displaystyle \underset{e}{\int }v{n}_2}}\right)}^2\right\}}\\ \le {\displaystyle \underset{i=1}{\overset{K}{\sum }}\left\{{\left({\displaystyle \underset{e\in {\epsilon }_h\left({\text{Γ}}_i\right)}{\sum }\left|n_1\right|\sqrt{\left|e\right|}{\left({\displaystyle \underset{e}{\int }{v}^2}\right)}^{\frac{1}{2}}}\right)}^2+{\left({\displaystyle \underset{e\in {\epsilon }_h\left({\text{Γ}}_i\right)}{\sum }\left|n_2\right|\sqrt{\left|e\right|}{\left({\displaystyle \underset{e}{\int }{v}^2}\right)}^{\frac{1}{2}}}\right)}^2\right\}}\\ \le \underset{i=1,\cdots ,K}{\max }\left|{\text{Γ}}_i\right|{\displaystyle \underset{e\in {\epsilon }_h\left(\Gamma \right)}{\sum }{\displaystyle \underset{e}{\int }{v}^2}}=\underset{i=1,\cdots ,K}{\max }\left|{\text{Γ}}_i\right|{\displaystyle \underset{e\in {\epsilon }_h\left(\Gamma \right)}{\sum }{\Vert v\Vert }_{L^2\left(e\right)}^2}\\ \le C{\Vert v\Vert }_{1,\text{Ω}}^2\le C{\left|v\right|}_{1,\text{Ω}}^2,\end{array}$

因此 $\sqrt{b\left(\cdot ,\cdot \right)}\lesssim {\left|\cdot \right|}_{1,\Omega }$。

相应于(2.2)的离散特征值问题为：求 ${\lambda }_h\in ℝ$ 和 $u_h\in V_h$ 满足

$\{\begin{array}{l}a\left(u_h,v_h\right)={\lambda }_{h}b\left(u_h,v_h\right),\forall v_h\in V_h,\hfill \\ b\left(u_h,u_h\right)=1.\hfill \end{array}$ (2.3)

(2.3)归结为一个广义矩阵特征值问题。由文献 [27] 可知，离散问题(2.3)具有2K个正的特征值：

$0<{\lambda }_{1,h}\le \cdots \le {\lambda }_{2K,h}$，每个特征值都对应一个特征函数 $u_j\in V$，使得 $\left\{u_{1,h},\cdots ,u_{2K,h}\right\}$ 是一个线性无关的集合。

相应于(2.2)和(2.3)的源问题和近似源问题分别如下：

求 $w\in V$，使得

$a\left(w,v\right)=b\left(f,v\right),\forall v\in V\text{.}$ (2.4)

求 $w_h\in V_h$，使得

$\begin{array}{r}\hfill a\left(w_h,v\right)=b\left(f,v\right),\forall v\in V_h\text{.}\end{array}$ (2.5)

由Neuman问题的标准先验估计(参见文献 [28] )可知，问题(2.4)的解 $\begin{array}{r}\hfill w\in H^{1+r}\left(\Omega \right)\end{array}$，其中 $r<\frac{\pi }{\theta }$， $\theta$ 为 $\Omega$ 的最大凹角，且有

$\begin{array}{r}\hfill {\Vert w\Vert }_{H^{1+r}\left(\Omega \right)/ℝ}\le C{\left|f\right|}_{1,\Omega }\text{.}\end{array}$ (2.6)

故对所有的 $r<\frac{\pi }{\theta }$，(2.2)的特征函数也满足 $u\in H^{1+r}\left(\Omega \right)$。注意到对于至少有一个交界面 ${\text{Γ}}_i$ 的多角形区域 $\Omega$，必有 $\theta >\pi$，因此 $\frac{1}{2}<\frac{\pi }{\theta }<1$。

由Lax-Milgram定理可知，(2.4)和(2.5)分别有唯一解，由此可定义下述有界线性算子 $T,T_h:V\to V$ ：

$\begin{array}{r}\hfill \{\begin{array}{l}f\in V\mapsto Tf\in V,\hfill \\ a\left(Tf,v\right)=b\left(f,v\right),\forall v\in V\hfill \end{array}\end{array}$

$\begin{array}{r}\hfill \{\begin{array}{l}f\in V\mapsto T_{h}f\in V_h\subset V,\hfill \\ a\left(T_{h}f,v_h\right)=b\left(f,v_h\right),\forall v_h\in V_h\hfill \end{array}\end{array}$

对任意的 $w\in V$，如下定义到 $V_h$ 的V椭圆投影算子 $P_h$ ：

$\begin{array}{r}\hfill P_{h}w\in V_h:a\left(P_{h}w-w,v_h\right)=0,\forall v_h\in V_h.\end{array}$

于是，对 $\forall f\in V$ 有

$\begin{array}{r}\hfill a\left(T_{h}f-P_h\left(Tf\right),v_h\right)=a\left(T_{h}f-Tf+Tf-P_h\left(Tf\right),v_h\right)=0,\forall v_h\in V_h\text{.}\end{array}$

故对 $\forall f\in V$，有 $T_{h}f=P_{h}Tf$，从而 $T_h=P_{h}T$。

由插值误差估计和(2.6)可推出

$\begin{array}{c}{\left|T_h-T\right|}_{1,\Omega }=\underset{g\in V}{\sup }\frac{{\left|\left(T_h-T\right)g\right|}_{1,\Omega }}{{\left|g\right|}_{1,\Omega }}=\underset{g\in V}{\sup }\frac{{\left|P_{h}Tg-Tg\right|}_{1,\Omega }}{{\left|g\right|}_{1,\Omega }}\\ \lesssim \underset{g\in V}{\sup }\frac{h^r{\left|Tg\right|}_{1+r,\Omega }}{{\left|g\right|}_{1,\Omega }}\lesssim \underset{g\in V}{\sup }\frac{h^r{\left|g\right|}_{1,\Omega }}{{\left|g\right|}_{1,\Omega }}=h^r\to 0\left(h\to 0\right),\end{array}$

从而得到下述引理。

引理2.1 ${\left|T-T_h\right|}_{1,\Omega }\to 0\left(h\to 0\right)$。

参考文献 [25] 我们有以下协调有限元逼近的误差估计。

引理2.2对于任意 $r<\frac{\pi }{\theta }$，存在正常数 $\kappa$，使得若 $h<\kappa$，则有

${\left|u_j-u_{j,h}\right|}_{1,\Omega }\lesssim h^r,$ (2.7)

$b{\left(u_j-u_{j,h},u_j-u_{j,h}\right)}^{\frac{1}{2}}\lesssim h^r{\left|u_j-u_{j,h}\right|}_{1,\Omega },$ (2.8)

$\left|{\lambda }_j-{\lambda }_{j,h}\right|\lesssim {\left|u_j-u_{j,h}\right|}_{1,\Omega }^2.$ (2.9)

由文献 [29] [30] 可知，(2.2)和(2.3)分别有如下等价的算子形式：

$Tu=\mu u,$

$T_h{u}_h={\mu }_h{u}_h,$

其中 $\mu =\frac{1}{\lambda }$， ${\mu }_h=\frac{1}{{\lambda }_h}$。在本文中 $\mu$ 和 ${\mu }_h$， $\lambda$ 和 ${\lambda }_h$ 都被称作特征值。

设 ${\mu }_k$ 的代数重数为q， ${\mu }_k={\mu }_{k+1}=\cdots ={\mu }_{k+q-1}$。令 $M\left({\mu }_k\right)$ 是T的所有对应于 ${\mu }_k$ 的特征函数张成的空间， $M_h\left({\mu }_k\right)$ 是 $T_h$ 的所有收敛到 ${\mu }_k$ 的特征值的特征空间的直和。令 $\hat{M}\left({\mu }_k\right)=\left\{v:v\in M\left({\mu }_k\right),{\left|v\right|}_{1,\Omega }=1\right\}$， ${\hat{M}}_h\left({\mu }_k\right)=\left\{v:v\in M_h\left({\mu }_k\right),{\left|v\right|}_{1,\Omega }=1\right\}$。我们也记 $M\left({\lambda }_k\right)=M\left({\mu }_k\right)$， $M_h\left({\lambda }_k\right)=M_h\left({\mu }_k\right)$， $\hat{M}\left({\lambda }_k\right)=\hat{M}\left({\mu }_k\right)$， ${\hat{M}}_h\left({\lambda }_k\right)={\hat{M}}_h\left({\mu }_k\right)$。

记

${\delta }_h\left({\lambda }_k\right)=\underset{w\in \hat{M}\left({\lambda }_k\right)}{\sup }\underset{v\in V_h}{\inf }{\left|w-v\right|}_{1,\Omega }.$

由(2.7)可知 ${\delta }_h\left({\lambda }_k\right)\lesssim h^r$。

引理2.3令 ${\lambda }_k$ 和 ${\lambda }_{k,h}$ 分别是(2.2)和(2.3)的第k个特征值，则

$\left|{\lambda }_{k,h}-{\lambda }_k\right|\lesssim {\delta }_h^2\left({\lambda }_k\right),$ (2.10)

对任意相应于 ${\lambda }_{k,h}$ 的特征函数 $u_{k,h}$ 且 ${\left|u_{k,h}\right|}_{1,\Omega }=1$，存在 $u_k\in \hat{M}\left({\lambda }_k\right)$ 使得

${\left|u_{k,h}-u_k\right|}_{1,\Omega }\lesssim {\delta }_h\left({\lambda }_k\right),$ (2.11)

对任意 $u_k\in \hat{M}\left({\lambda }_k\right)$，存在 $u_h\in M_h\left({\lambda }_k\right)$ 使得

${\left|u_h-u_k\right|}_{1,\Omega }\lesssim {\delta }_h\left({\lambda }_k\right).$ (2.12)

引理2.4令 $\left(\lambda ,u\right)$ 是(2.2)的特征对，则对任意的 $w\in V$， $b\left(w,w\right)\ne 0$，Rayleigh商 $\frac{a\left(w,w\right)}{b\left(w,w\right)}$ 满足

$\begin{array}{r}\hfill \frac{a\left(w,w\right)}{b\left(w,w\right)}-\lambda =\frac{{\left|w-u\right|}_{1,\Omega }}{b\left(w,w\right)}-\lambda \frac{b\left(w-u,w-u\right)}{b\left(w,w\right)}\end{array}.$ (2.13)

证明：详见文献 [30] 中引理 9.1。

下面的引理(参见文献 [13] 中定理3.2)是分析多网格离散方案的基本工具。

引理2.5令 $\left({\mu }_0,u_0\right)$ 是 $T_h$ 的一个近似特征对， $u_0\in V_h$ 且 ${\left|u_0\right|}_{1,\Omega }=1$。假设 $\underset{k\le j\le k+q+1}{\max }\left|\frac{{\mu }_{j,h}-{\mu }_{k,h}}{{\mu }_0-{\mu }_{k,h}}\right|\le 1/2$， $dist\left(u_0,M_h\left({\mu }_k\right)\right)\le 1/2$， $\left|{\mu }_0-{\mu }_{j,h}\right|\ge \vartheta /2\left(j\ne k,k+1,\cdots ,k+q-1\right)$， $u\in V_h$ 与 $u_k^h\in V_h$ 满足

$\left({\mu }_0-T_h\right)u=u_0,u_k^h=\frac{u}{{\left|u\right|}_{1,\Omega }},$

则

$dist\left(u_k^h,{\hat{M}}_h\left({\mu }_k\right)\right)\le \frac{16}{\vartheta }\left|{\mu }_0-{\mu }_{k,h}\right|dist\left(u_0,M_h\left({\mu }_k\right)\right),$ (2.14)

其中 $\vartheta =\underset{{\mu }_j\ne {\mu }_k}{\min }\left|{\mu }_j-{\mu }_k\right|$ 是特征值 ${\mu }_k$ 的分隔常数。

3. 多网格方案

本节中，我们结合有限元方法和Rayleigh商移位反迭代法建立流固振动Laplace模型问题的多网格离散方案。令 ${\left\{{\pi }_{h_i}\right\}}_{i=1}^l$ 是一族正规网格， $h_{i-1}\gg h_i$，且令 $V_{h_i}$ 是定义在 ${\pi }_{h_i}$ 上的协调有限元空间，令 ${\pi }_H={\pi }_{h_1}$， $V_H=V_{h_1}$， ${\pi }_h={\pi }_{h_l}$， $V_h=V_{h_l}$。

方案1多网格方案

步骤1在初始粗网格 ${\pi }_H$ 上解(2.3)：求 ${\lambda }_{k,H}\in ℝ$， $u_{k,H}\in V_H$ 使得 ${\left|u_{k,H}\right|}_{1,\Omega }=1$，且

$a\left(u_{k,H},v\right)={\lambda }_{k,H}b\left(u_{k,H},v\right),\forall v\in V_H.$

步骤2执行赋值命令： $u_k^{h_1}\Leftarrow u_{k,H}$， ${\lambda }_k^{h_1}\Leftarrow {\lambda }_{k,H}$， $i\Leftarrow 2$。

步骤3在 ${\pi }_{h_i}$ 上解线性方程组：求 $\tilde{u}\in V_{h_i}$ 使得

$a\left(\tilde{u},v\right)-{\lambda }_k^{h_{i-1}}b\left(\tilde{u},v\right)=b\left(u_k^{h_{i-1}},v\right),\forall v\in V_{h_i},$

且取 $u_k^{h_i}=\tilde{u}/{\left|\tilde{u}\right|}_{1,\Omega }$。

步骤4计算Rayleigh商

${\lambda }_k^{h_i}=\frac{a\left(u_k^{h_i},u_k^{h_i}\right)}{b\left(u_k^{h_i},u_k^{h_i}\right)}.$

步骤5如果 $i=l$，则输出 $\left({\lambda }_k^{h_l},u_k^{h_l}\right)$，即输出 $\left({\lambda }_k^h,u_k^h\right)$，停止。否则， $i\Leftarrow i+1$，返回步骤3。

下面分析由方案1求得的近似解的误差。

定理3.1.令 $\left({\lambda }_k^{h_l},u_k^{h_l}\right)$ 是方案1所求得的近似特征对。假设 $\left({\lambda }_k^{h_{l-1}},u_k^{h_{l-1}}\right)$ 逼近(2.2)的一个特征对 $\left({\lambda }_k,\bar{u}\right)$， $\bar{u}\in \hat{M}\left({\lambda }_k\right)$，且 ${\lambda }_k^{h_{l-1}}-{\lambda }_k$ 相较于 ${\lambda }_{k,h_l}-{\lambda }_k$ 是一个低阶的无穷小量，则存在 $u_k\in M\left({\lambda }_k\right)$ 使得

${\left|u_k^{h_l}-u_k\right|}_{1,\Omega }\lesssim \frac{48}{\vartheta }{\lambda }_k\left|{\lambda }_k^{h_{l-1}}-{\lambda }_k\right|dist\left(u_k^{h_{l-1}},\hat{M}\left({\lambda }_k\right)\right)+3q{\delta }_{h_l}\left({\lambda }_k\right),$ (3.1)

$\left|{\lambda }_k^{h_l}-{\lambda }_k\right|\lesssim 2{\lambda }_k\left(1+{\lambda }_k\right){\left|u_k^{h_l}-u_k\right|}_{1,\Omega }^2,l\ge 2.$ (3.2)

证明：我们利用引理2.5来完成证明。首先验证引理2.5的条件是满足的。

选取 ${\mu }_0=1/{\lambda }_k^{h_{l-1}}$， $u_0=\frac{{\lambda }_k^{h_{l-1}}{T}_{h_l}{u}_k^{h_{l-1}}}{{\left|{\lambda }_k^{h_{l-1}}{T}_{h_l}{u}_k^{h_{l-1}}\right|}_{1,\Omega }}$。由引理2.1可知 ${\left|T_{h_l}-T\right|}_{1,\Omega }\to 0\left(h_l\to 0\right)$，则有

${\left|T_{h_l}f\right|}_{1,\Omega }\lesssim {\left|f\right|}_{1,\Omega },\forall f\in V.$

因此，由定理假设可推出

$\begin{array}{l}{\left|{\lambda }_k^{h_{l-1}}{T}_{h_l}{u}_k^{h_{l-1}}-\bar{u}\right|}_{1,\Omega }\\ ={\left|{\lambda }_k^{h_{l-1}}{T}_{h_l}{u}_k^{h_{l-1}}-{\lambda }_k^{h_{l-1}}{T}_{h_l}\bar{u}+{\lambda }_k^{h_{l-1}}{T}_{h_l}\bar{u}-{\lambda }_k{T}_{h_l}\bar{u}+{\lambda }_k{T}_{h_l}\bar{u}-{\lambda }_{k}T\bar{u}\right|}_{1,\Omega }\\ \lesssim {\lambda }_k^{h_{l-1}}{\left|u_k^{h_{l-1}}-\bar{u}\right|}_{1,\Omega }+\left|{\lambda }_k^{h_{l-1}}-{\lambda }_k\right|{\left|\bar{u}\right|}_{1,\Omega }+{\lambda }_k{\left|\left(T_{h_l}-T\right)\bar{u}\right|}_{1,\Omega }.\end{array}$

注意到对任一赋范空间中的任意非零函数 $u,v$，成立下述不等式

$\Vert \frac{u}{\Vert u\Vert }-\frac{v}{\Vert v\Vert }\Vert \le 2\frac{\Vert u-v\Vert }{\Vert u\Vert },\Vert \frac{u}{\Vert u\Vert }-\frac{v}{\Vert v\Vert }\Vert \le 2\frac{\Vert u-v\Vert }{\Vert v\Vert }.$ (3.3)

因此有

$\begin{array}{c}dist\left(u_0,\hat{M}\left({\lambda }_k\right)\right)\le {\left|u_0-\bar{u}\right|}_{1,\Omega }\le 2{\left|{\lambda }_k^{h_{l-1}}{T}_{h_l}{u}_k^{h_{l-1}}-\bar{u}\right|}_{1,\Omega }\\ \lesssim 2\left\{{\lambda }_k^{h_{l-1}}{\left|u_k^{h_{l-1}}-\bar{u}\right|}_{1,\Omega }+\left|{\lambda }_k^{h_{l-1}}-{\lambda }_k\right|{\left|\bar{u}\right|}_{1,\Omega }+{\lambda }_k{\left|\left(T_{h_l}-T\right)\bar{u}\right|}_{1,\Omega }\right\}.\end{array}$ (3.4)

由三角不等式和(2.12)可得

$dist\left(u_0,M_{h_l}\left({\lambda }_k\right)\right)\lesssim dist\left(u_0,\hat{M}\left({\lambda }_k\right)\right)+{\delta }_{h_l}\left({\lambda }_k\right).$ (3.5)

由(2.10)可知 ${\lambda }_{k,h_l}\to {\lambda }_k\left(h_l\to 0\right)$，于是由定理假设有

$\left|{\mu }_0-{\mu }_{k,h_l}\right|=\left|\frac{{\lambda }_k^{h_{l-1}}-{\lambda }_k+{\lambda }_k-{\lambda }_{k,h_l}}{{\lambda }_{k,h_l}{\lambda }_k^{h_l}}\right|=O\left(\left|{\lambda }_k^{h_{l-1}}-{\lambda }_k\right|\right).$ (3.6)

由(2.5)可知方案1中的步骤3等价于

$a\left(\tilde{u},v\right)-{\lambda }_k^{h_{l-1}}a\left(T_{h_l}\tilde{u},v\right)=a\left(T_{h_l}{u}_k^{h_{l-1}},v\right),\forall v\in V_{h_l},u_k^{h_l}=\frac{\tilde{u}}{{\left|\tilde{u}\right|}_{1,\Omega }},$

即

$\left(\frac{1}{{\lambda }_k^{h_{l-1}}}-T_{h_l}\right)\tilde{u}=\frac{1}{{\lambda }_k^{h_{l-1}}}{T}_{h_l}{u}_k^{h_{l-1}},u_k^{h_l}=\frac{\tilde{u}}{{\left|\tilde{u}\right|}_{1,\Omega }}\text{.}$ (3.7)

注意到 $\frac{1}{{\lambda }_k^{h_{l-1}}}{T}_{h_l}{u}_k^{h_{l-1}}={\left|\frac{1}{{\lambda }_k^{h_{l-1}}}{T}_{h_l}{u}_k^{h_{l-1}}\right|}_{1,\Omega }{u}_0$ 与 $u_0$ 仅仅相差一个常数倍，从而方案1步骤3等价于

$\left(\frac{1}{{\lambda }_k^{h_{l-1}}}-T_{h_l}\right)\tilde{u}=u_0,u_k^{h_l}=\frac{\tilde{u}}{{\left|\tilde{u}\right|}_{1,\Omega }}.$ (3.8)

事实上，由(3.7)，(3.8)求得的 $u_k^{h_l}$ 显然是相等的。当 $h_{l-1}$ 足够小时，注意到 $h_l\ll h_{l-1}$，由(3.4)和(3.5)得

$dist\left(u_0,M_{h_l}\left({\lambda }_k\right)\right)\le \frac{1}{2}.$

因为 ${\lambda }_k={\lambda }_{k+1}=\cdots ={\lambda }_{k+q-1}$，故有

$\left|{\mu }_{j,h_l}-{\mu }_{k,h_l}\right|=\left|\frac{{\lambda }_{k,h_l}-{\lambda }_{j,h_l}}{{\lambda }_{k,h_l}{\lambda }_{j,h_l}}\right|=\left|\frac{{\lambda }_{k,h_l}-{\lambda }_k+{\lambda }_j-{\lambda }_{j,h_l}}{{\lambda }_{k,h_l}{\lambda }_{j,h_l}}\right|,$ (3.9)

且由(2.10)可得

${\lambda }_{j,h_l}-{\lambda }_k\lesssim {\delta }_{h_l}^2\left({\lambda }_k\right)\left(j=k,k+1,\cdots ,k+q-1\right).$

结合(3.9)和(3.6)，且注意到(3.9)的右端项相较于 ${\lambda }_k^{h_l}-{\lambda }_k$ 是高阶无穷小量，于是

$\underset{k\le j\le k+q-1}{\max }\left|\frac{{\mu }_{j,h_l}-{\mu }_{k,h_l}}{{\mu }_0-{\mu }_{j,h_l}}\right|\le \frac{1}{2}.$

因为 $\vartheta$ 是分隔常数， $h_{l-1}$ 足够小，且 $h_l\ll h_{l-1}$，故成立

$\left|{\mu }_0-{\mu }_{j,h_l}\right|\ge \frac{\vartheta }{2},j\ne k,k+1,\cdots ,k+q-1.$

综上所述，可知引理2.5的条件是成立的。

下面证明(3.1)和(3.2)成立。

将(3.5)和(3.6)代入(2.14)，得到

$dist\left(u_k^{h_l},{\hat{M}}_{h_l}\left({\mu }_k\right)\right)\lesssim \frac{16}{\vartheta }\left|{\lambda }_k^{h_{l-1}}-{\lambda }_k\right|\left(dist\left(u_0,\hat{M}\left({\lambda }_k\right)\right)+{\delta }_{h_l}\left({\lambda }_k\right)\right).$ (3.10)

令特征向量 ${\left\{u_{j,h_l}\right\}}_{j=k}^{k+q-1}$ 是 $M_{h_l}\left({\lambda }_k\right)$ 的在内积 $a\left(\cdot ,\cdot \right)$ 意义下的一组正交基，并注意到

$dist\left(u_k^{h_l},M_{h_l}\left({\lambda }_k\right)\right)={\left|u_k^{h_l}-{\displaystyle \underset{j=k}{\overset{k+q-1}{\sum }}a\left(u_k^{h_l},u_{j,h_l}\right)u_{j,h_l}}\right|}_{1,\Omega }.$

令

$u^{\text{*}}={\displaystyle \underset{j=k}{\overset{k+q-1}{\sum }}a}\left(u_k^{h_l},u_{j,h_l}\right)u_{j,h_l},$

则由(3.10)可推得

${\left|u_k^{h_l}-u^{\text{*}}\right|}_{1,\Omega }\lesssim \frac{16}{\vartheta }\left|{\lambda }_k^{h_{l-1}}-{\lambda }_k\right|\left(dist\left(u_0,\hat{M}\left({\lambda }_k\right)\right)+{\delta }_{h_l}\left({\lambda }_k\right)\right).$ (3.11)

由引理2.3知，存在 ${\left\{u_j^0\right\}}_{j=k}^{k+q-1}\subset \hat{M}\left({\lambda }_k\right)$ 使得 $u_{j,h_l}-u_j^0$ 满足(2.11)。令

$u_k={\displaystyle \underset{j=k}{\overset{k+q-1}{\sum }}a}\left(u_k^{h_l},u_{j,h_l}\right)u_j^0,$

则 $u_k\in M\left({\lambda }_k\right)$。由(2.11)可推得

$\begin{array}{c}{\left|u^{*}-u_k\right|}_{1,\Omega }={\left|{\displaystyle \underset{j=k}{\overset{k+q-1}{\sum }}a}\left(u_k^{h_l},u_{j,h_l}\right)\left(u_{j,h_l}-u_j^0\right)\right|}_{1,\Omega }\lesssim {\left({\displaystyle \underset{j=k}{\overset{k+q-1}{\sum }}{\left|u_{j,h_l}-u_j^0\right|}_{1,\Omega }^2}\right)}^{\frac{1}{2}}\\ \lesssim {\displaystyle \underset{j=k}{\overset{k+q-1}{\sum }}{\delta }_{h_l}\left({\lambda }_j\right)}\lesssim q{\delta }_{h_l}\left({\lambda }_k\right).\end{array}$

结合(3.11)与上述不等式可得

$\begin{array}{c}{\left|u_k^{h_l}-u_k\right|}_{1,\Omega }\le {\left|u_k^{h_l}-u^{\text{*}}\right|}_{1,\Omega }+{\left|u^{\text{*}}-u_k\right|}_{1,\Omega }\\ \lesssim \frac{16}{\vartheta }\left|{\lambda }_k^{h_{l-1}}-{\lambda }_k\right|\left(dist\left(u_0,\hat{M}\left({\lambda }_k\right)\right)+{\delta }_{h_l}\left({\lambda }_k\right)\right)+q{\delta }_{h_l}\left({\lambda }_k\right)\\ \lesssim \frac{16}{\vartheta }\left|{\lambda }_k^{h_{l-1}}-{\lambda }_k\right|dist\left(u_0,\hat{M}\left({\lambda }_k\right)\right)+2q{\delta }_{h_l}\left({\lambda }_k\right).\end{array}$ (3.12)

显然存在 ${u^{\prime }}_k\in \hat{M}\left({\lambda }_k\right)$ 使得

${\left|u_k^{h_{l-1}}-{u^{\prime }}_k\right|}_{1,\Omega }=dist\left(u_k^{h_{l-1}},\hat{M}\left({\lambda }_k\right)\right),$

${\lambda }_{k}T{u^{\prime }}_k={u^{\prime }}_k\in \hat{M}\left({\lambda }_k\right)$，以及 ${\left|\left(T_{h_l}-T\right){u^{\prime }}_k\right|}_{1,\Omega }\le \frac{1}{{\lambda }_k}{\delta }_{h_l}\left({\lambda }_k\right)$。

因为 $u_0={\lambda }_k^{h_{l-1}}{T}_{h_l}{u}_k^{h_{l-1}}/{\left|{\lambda }_k^{h_{l-1}}{T}_{h_l}{u}_k^{h_{l-1}}\right|}_{1,\Omega }$，故由(2.13)和(3.3)推出

$\begin{array}{c}dist\left(u_0,\hat{M}\left({\lambda }_k\right)\right)\lesssim 2dist\left({\lambda }_k^{h_{l-1}}{T}_{h_l}{u}_k^{h_{l-1}},\hat{M}\left({\lambda }_k\right)\right)\le 2{\left|{\lambda }_k^{h_{l-1}}{T}_{h_l}{u}_k^{h_{l-1}}-{\lambda }_{k}T{u^{\prime }}_k\right|}_{1,\Omega }\\ \le 2{\left|\left({\lambda }_k^{h_{l-1}}-{\lambda }_k\right)T_{h_l}{u}_k^{h_{l-1}}\right|}_{1,\Omega }+2{\left|{\lambda }_k{T}_{h_l}\left(u_k^{h_{l-1}}-{u^{\prime }}_k\right)\right|}_{1,\Omega }+2{\left|{\lambda }_k\left(T_{h_l}-T\right){u^{\prime }}_k\right|}_{1,\Omega }\\ \lesssim 2\left|{\lambda }_k^{h_{l-1}}-{\lambda }_k\right|+2{\lambda }_{k}dist\left(u_k^{h_{l-1}},\hat{M}\left({\lambda }_k\right)\right)+2{\delta }_{h_l}\left({\lambda }_k\right)\\ \lesssim 3{\lambda }_{k}dist\left(u_k^{h_{l-1}},\hat{M}\left({\lambda }_k\right)\right)+2{\delta }_{h_l}\left({\lambda }_k\right).\end{array}$ (3.13)

将(3.13)代入(3.12)，便得到(3.1)。

仍记 $\frac{u_k}{{\left|u_k\right|}_{1,\Omega }}$ 为 $u_k$，则联系(3.1)与(3.3)可知 $u_k^{h_l}$ 在 ${\left|\cdot \right|}_{1,\text{Ω}}$ 意义下收敛到 $u_k$。由于 $\sqrt{b\left(\cdot ,\cdot \right)}\lesssim {\left|\cdot \right|}_{1,\text{Ω}}$，故 $b\left(u_k^{h_l}-u_k,u_k^{h_l}-u_k\right)\to 0$，于是 $\frac{1}{b\left(u_k^{h_l},u_k^{h_l}\right)}\to \frac{1}{b\left(u_k,u_k\right)}={\lambda }_k$。因此，当 $h_l$ 足够小时，成立 $\frac{1}{b\left(u_k^{h_l},u_k^{h_l}\right)}\le 2{\lambda }_k$。在(2.13)中取 $w=u_k^{h_l}$，得到

$\begin{array}{c}\left|{\lambda }_k^{h_l}-{\lambda }_k\right|\le \frac{{\left|u_k^{h_l}-u_k\right|}_{1,\Omega }^2}{b\left(u_k^{h_l},u_k^{h_l}\right)}+{\lambda }_k\frac{b\left(u_k^{h_l}-u_k,u_k^{h_l}-u_k\right)}{b\left(u_k^{h_l},u_k^{h_l}\right)}\\ \le 2{\lambda }_k{\left|u_k^{h_l}-u_k\right|}_{1,\Omega }^2+2{\lambda }_k^{2}b\left(u_k^{h_l}-u_k,u_k^{h_l}-u_k\right)\\ \lesssim 2{\lambda }_k\left(1+{\lambda }_k\right){\left|u_k^{h_l}-u_k\right|}_{1,\Omega }^2.\end{array}$

证毕。

在方案1中取 $l=2$，并记 $H=h_1$， $h=h_2$，我们立即得到下面的二网格方案。

方案2二网格方案

步骤1 在 ${\pi }_H$ 上解(2.3)：求 ${\lambda }_{k,H}\in ℝ$， $u_{k,H}\in V_H$ 使得 ${\left|u_{k,H}\right|}_{1,\Omega }=1$，且

$a\left(u_{k,H},v\right)={\lambda }_{k,H}b\left(u_{k,H},v\right),\forall v\in V_H.$

步骤2在 ${\pi }_h$ 上解线性方程组：求 $\tilde{u}\in V_h$ 使得

$\begin{array}{c}a\left(\tilde{u},v\right)-{\lambda }_{k,H}b\left(\tilde{u},v\right)=b\left(u_{k,H},v\right),\forall v\in V_h,\end{array}$

且取 $u_k^h=\tilde{u}/{\left|\tilde{u}\right|}_{1,\Omega }$。

步骤3计算Rayleigh商

${\lambda }_k^h=\frac{a\left(u_k^h,u_k^h\right)}{b\left(u_k^h,u_k^h\right)}\text{.}$

由定理3.1可以得到二网格离散方案的误差估计如下：

${\left|u_k^h-u_k\right|}_{1,\Omega }\lesssim \frac{48}{\vartheta }{\lambda }_k\left|{\lambda }_{k,H}-{\lambda }_k\right|dist\left(u_{k,H},\hat{M}\left({\lambda }_k\right)\right)+3q{\delta }_h\left({\lambda }_k\right),$

$\left|{\lambda }_k^h-{\lambda }_k\right|\lesssim 2{\lambda }_k\left(1+{\lambda }_k\right){\left|u_k^h-u_k\right|}_{1,\Omega }^2\text{.}$

注意到

$\left|{\lambda }_{k,H}-{\lambda }_k\right|\lesssim {\delta }_H^2\left({\lambda }_k\right)\lesssim H^{2r},$

$dist\left(u_{k,H},\hat{M}\left({\lambda }_k\right)\right)\lesssim {\delta }_H\left({\lambda }_k\right)\lesssim H^r,$

因此

${\left|u_k^h-u_k\right|}_{1,\Omega }\lesssim H^{3r}+h^r,$ (3.14)

$\left|{\lambda }_k^h-{\lambda }_k\right|\lesssim H^{6r}+h^{2r}.$ (3.15)

(3.15)表明，当 $H=O\left(\sqrt[3]{h}\right)$ 时，由方案2求得的近似特征值可以达到渐进最优收敛阶。

4. 数值实验

本节将呈现一些数值算例来展示本文建立的多网格离散方案的有效性。数值实验是在具有1.8 GHZ CPU和4 G RAM的DellInspiron 5547 PC上，在MATLAB 2015a环境中借助软件包IFEM (参见文献 [31] )进行编程计算的。我们利用方案2，采取线性协调元完成数值算例，在实验表格中采用下述符号：

H：粗网格 ${\pi }_H$ 的直径。

h：细网格 ${\pi }_h$ 的直径。

${\lambda }_{j,H}$ ：在网格 ${\pi }_H$ 上利用eigs命令直接求解得到的(2.3)第j个特征值。

${\lambda }_j^h$ ：利用方案2求得的第j个特征值。

$tim{e}_h\left(秒\right)$ ：直接在细网格 ${\pi }_h$ 上求解特征值问题所用的cpu时间。

$tim{e}^h\left(秒\right)$ ：用方案2计算，从程序开始运行至当前结果出现的cpu时间。

−：表示由于计算机内存限制无法计算。

算例1. 在一个边长为8的四边形方形腔，其中心有一个边长为 $2\sqrt{2}$ 的菱形管的区域上考虑问题(2.1)，如图2 (左) 所示。计算结果列在表1中，同时还列出了直接求解特征值问题所得的结果。

文献 [25] 中给出了问题(2.1)在此区域上最小特征值的参考值： ${\lambda }_1\approx 0.07896$，且其代数重数为2。由表1可以看到，与直接在细网格上求解相比较，利用方案2可以用较少的时间求得精度较高的近似特征值。

Figure 2. Square cavity with a rhomboidal tube (left) and Square cavity with two square tubes (right)

图2. 内部有一个菱形管的方形腔(左)，内部有两个方形管的方形腔(右)

Table 1. The approximate eigenvalue on the square cavity with a rhomboidal tube

表1. 内部有一个菱形管的方形腔上的近似特征值

算例2在一个边长为8的四边形方形腔，其内部有两个边长为2的方形管的区域上考虑问题(2.1)，如图2(右)所示。计算结果列在表2中，同时还列出了直接求解特征值问题所得的结果。由表2可以看到，利用方案2可以用较少时间求得与直接在细网格上计算所得特征值精度相同的近似解，而且直接计算无法进行时，利用方案2仍然可以计算。

Table 2. The approximate eigenvalueson the square cavity with two square tubes

表2. 内部有两个方形管的方形腔上的近似特征值

基金项目

国家自然科学基金项目(批准号：11761022)。