谱元法求解Helmholtz方程透射特征值问题*

摘要：　研究了Helmholtz方程透射特征值问题，提出一种Chebyshev谱元法求解，该方法兼具了有限元法处理边界及区域的灵活性和谱方法的快速收敛特性．运用加权余量原理，得到了Chebyshev谱元法用于透射特征值问题的基本理论以及数学公式，将原问题转化为二次特征值问题．最后通过数值实验算例验证了Chebyshev谱元法的有效性．

关　键　词：　透射特征值问题；　二次特征值问题；　谱元法；　Chebyshev基函数

引　　言

透射特征问题来源于非均匀的逆散射理论[1]，在近几年的发展中已经成为逆散射理论的一个重要组成部分[2]．特征值能够通过远场数据分析得到并且可以用来分析散射对象的材料属性[3]，在很多领域有广泛的应用．目前，对于透射特征值问题的数值模拟有比较多的方法，如有限元法、边界元法[4] ,Colton和Sun等对透射特征值问题提出了不同的离散方法[5-8]：混合有限元法、二次特征值法和连续有限元法等．有限元法能够很好地处理边界和区域复杂结构，有很好的适应性，但是随着对象波数的增大，需要更多的网络节点数目，这样会使得计算量增大，耗费大量计算资源．对于这类问题，谱方法是一种比较好的处理方法[9]．它将插值基函数取为无限可微全局函数，具有较好的收敛性．理论上谱方法可以在最短波长上取最少的网格点得到所需的精度[10],但是谱方法对于复杂区域上的问题求解效果不好，区域适应性不高．谱元法是由Patera[11]在结合谱方法的精度和有限元的思想的基础上提出来的一种新方法，它兼具了有限元处理边界和区域的灵活性和谱方法的快速收敛特性．由于在单元上采用无穷光滑的插值基函数，谱元法的精度可以根据需求和计算条件来调节[12]，并且具有和有限元同样优势的几何区域适应性，所以该方法非常适应求解透射特征值问题．

本文结构如下：第1节通过变分原理得到原问题的弱形式，并对二维有界单联通区域D进行了单元划分，构造Chebyshev谱元逼近函数[13]，得到总刚度矩阵和总荷载矩阵；第2节给出多组数值实验，通过对比其他方法，验证了谱元法的有效性；第3节对文章进行了小结和展望．

1　二维透射特征值问题数值方法

1.1　数值模型

考虑声波在一个有界单连通的各向同性介质D∈R2上的Helmholtz方程透射特征值问题，即求解

满足方程：

其中γ为散射体边界上的单位法向量，n(x)表示折射率，非零常数k为传输特征值，相应的w,v称为特征函数，称(k,w,v)为特征对．

1.2　标准单元内函数Chebyshev谱逼近

将有界单连通区域D分解为若干个互不重叠的矩形区域，记单元总数为

考察第i个单元Di,并设

将此单元Di化成标准正方形求解单元ei，则

用谱方法进行方程离散时，通常使用Legendre正交多项式或者Chebyshev正交多项式构造全局基函数[15]．本文选择Chebyshev正交多项式配置插值点，在单元ei内,在x方向选取

个点,在y方向选取

个点作为u(x,y)的插值点,此处用Chebyshev多项式的极值点作为插值点：

1.3　单元刚度矩阵的生成

将插值函数ui(ξ,η),vi(ξ,μ)代入式(6)中，并运用

的任意性，可得到

表1　标准单元区域上几种矩阵定义

Table 1　Several matrix definitions in the standard cell area

1.4　总体矩阵的生成

将各单元上矩阵合成为总体矩阵，最终可以得到以下一般二次特征值问题：

2　数值分析

为了验证Chebyshev谱元法求解透射特征值问题的有效性．本文中考虑3种数学模型： ① 以原点为中心的边长为2的正方形； ② 以原点为中心的环型，外边界为边长为2的正方形，内边界为边长为1的正方形，环宽度为1； ③ L型，长为3，宽为3，L型区域宽度为1．区域划分时将整个区域分解为等规模的若干个矩形，取网格大小h=1/12，在单元矩形上利用Chebyshev多项式的极值点作为插值点，本文均选择3阶Chebyshev多项式进行离散，如图1所示．

图1　3次Chebyshev多项式配置插值点示意图
Fig. 1　Collocation points of the 3rd-order Chebyshev interpolation

本文考虑n(x)取常数和非常数两种情况，在3个模型区域上得到4个最小的正实透射特征值，结果如表2、3所示．

表2　n(x)=16时不同区域上所对应的4个最小的正实透射特征值

Table 2　The first real transmission eigenvalues of some different regions.The index of refractionn(x)=16

表3　n(x)=x1-x2时不同区域上所对应的4个最小的正实透射特征值

Table 3　The first real transmission eigenvalues of some different regions. The index ofrefractionn(x)=x1-x2

图2　正方形区域上总体矩阵(nz=967) 图3　环形区域上总体矩阵(nz=2 862) 图4　L形区域上总体矩阵(nz=1 262)
Fig. 2　The global mass matrix in the square area (nz=967)ring area (nz=2 862) Fig. 3　The global mass matrix in the Fig. 4　The global mass matrix in the L area (nz=1 262)

图5　正方形模型上谱元法与有限元法结果对比
Fig. 5　Comparison between the spectral element method and the finite element method in the square area

图6　环型模型上谱元法和有限元法结果对比图7　L型模型上谱元法与有限元法结果对比
Fig. 6　Comparison between the spectral element method the finite element method in the ring area Fig. 7　Comparison between the spectral element method and and the finite element method in the L area

图8　在相同基函数阶数下谱元法与有限元运行时间对比
Fig. 8　The time consumptions of the spectral method and the finite element method in the same order of Chebyshev primary function

图2～图4分别给出了3种区域上总体矩阵结构．需要说明的是，由于谱元法所生成的总体矩阵是对称的、稀疏的、正定的、带状的，因此采用基于MPI的并行运算[16]，对于总体矩阵采用分布式存储．该方法可以较大地节约运行内存，从而提高执行效率．本文所有数值实验都是在MATLAB R2014a,Windows 10，系统安装内存2.00 GB下实现，本文将Chebyshev谱元法与传统有限元法进行对比，通过误差分析得到谱元法在计算精度上相对于有限元法较高，如图5～7所示．图8给出了谱元法和有限元法执行时间对比．

通过以上数据以及对比图像反映，Chebyshev谱元法求解二维区域透射特征值问题是合理有效的，对区域适应性也比较强．在相同精度要求下，谱元法通过较稀疏的单元划分可减少计算时间和运行内存．

3　总结和发展

本文利用Chebyshev谱元法对二维区域上透射特征值问题进行处理，首先通过变分法将原问题变为二次特征值问题．将总体区域划分为若干规则矩形，单元上采用三阶精度的Chebyshev谱函数进行逼近，谱元法通过在每一个单元上应用谱展开从而结合了谱方法和有限元的优点，将每个单元上得到的矩阵集成为总体矩阵．最后再选取3个特定模型进行数值求解，通过实验结果验证了Chebyshev谱元法求解透射特征值问题的有效性．

今后的工作可以从以下几个方面展开：加密网格，并采用更高阶Chebyshev谱逼近；针对形状更为复杂的散射体进行计算分析；三维区域上Chebyshev谱元法推广．

[1]　 COLTOND, MONK P. The inverse scattering problem for time-hormonic acoustic waves in an inhomogeneous medium[J]. The Quarterly Journal of Mechanics and Applied Mathematics, 1988, 26(1): 323-350.

[2]　 COLTON D, PAIVARINTA L, SYLVESTER J. The interior transmission problem[J]. Inverse Problems & Imaging, 2017, 1(1): 13-28.

[3]　 CAKONI F, COLTON D, HADDARH. On the determination of Dirichlet or transmission eigenvalues from far field data[J]. Comptes Rendus Mathematique, 2012, 348(7): 379-383.

[4]　 CAKONI F, KRESS R. A boundary integral equation method for the transmission eigenvalue problem[J]. Applicable Analysis, 2016, 96(1): 23-38.

[5]　 COLTON D, MONK P, SUN J. Analytical and computational methods for transmission eigenvalues[J]. Inverse Problems, 2010, 26(26): 045011.

[6]　 JI X, SUN J, TURNER T. Algorithm 922: a mixed finite element method for Helmholtz transmission eigenvalues[J]. ACM Transactions on Mathematical Software, 2012, 38(4): 1-8.

[7]　 JI X, SUN J, XIE H. A multigrid method for Helmholtz transmission eigenvalue problems[J]. Journal of Scientific Computing, 2014, 60(2): 276-294.

[8]　 SUN J. Iterative methods for transmission eigenvalues[J]. Society for Industrial and Applied Mathematics, 2011, 49(49): 1860-1874.

[9]　 AN J, SHEN J. Spectral approximation to a transmission eigenvalue problem and its applications to an inverse problem[J]. Computers & Mathematics With Applications, 2015, 69(10): 1132-1143.

[10]　ORSZAG S A. Spectral methods for problems in complex geometries[J]. Journal of Computational Physics, 1980, 37(1): 70-92.

[11]　PATERA A T. A spectral element method for fluid dynamics: laminar flow in a channel expansion[J]. Journal of Computational Physics, 1984, 54(3): 468-488.

[12]　容志建, 许传炬. 基于张量乘积的快速谱元算法[J]. 数学研究, 2008, 41(3): 264-271.(RONG Zhijian, XU Chuanju. Tensor product based fast spectral element solves[J]. Journal of Mathematical Study, 2008, 41(3): 264-271.(in Chinese))

[13]　林伟军. 弹性波传播模拟的Chebyshev谱元法[J]. 声学学报, 2007,32(6): 525-533.(LIN Weijun. A Chebyshev spectral element method for elastic wave modeling[J]. Acta Acustica, 2007, 32(6): 525-533.(in Chinese))

[14]　周欣, 李铁香. Helmholtz方程透射特征值问题的数值算法[J]. 应用数学进展, 2016, 5(4): 683-694.(ZHOU Xin, LI Tiexiang. Numerical solution of transmission eigenvalue problems of Helmholtz equation[J]. Advances in Applied Mathematics, 2016, 5(4): 683-694. (in Chinese))

[15]　朱晓钢, 聂玉峰. 变系数分数阶对流扩散方程的一种算子矩阵方法[J]. 应用数学和力学, 2018, 39(1): 104-112.(ZHU Xiaogang, NIE Yufeng. An operational matrix method for fractional advection-diffusion equations with variable coefficients[J]. Applied Mathematics and Mechanics, 2018, 39(1): 104-112.(in Chinese))

[16]　朱昌允, 秦国良, 徐忠. Chebyshev谱元方法结合并行算法求解三维区域的Helmholtz方程[J]. 应用力学学报, 2012, 29(3): 247-251.(ZHU Changyun, QIN Guoliang, XU Zhong. Parallel Chebyshev spectral element method for Helmholtz equation in 3D domain[J]. Chinese Journal of Applied Mechanics, 2012, 29(3): 247-251.(in Chinese))

A Spectral Element Method for Transmission Eigenvalue Problems of the Helmholtz Equation

Abstract: A Chebyshev spectral element method for the transmission eigenvalue problems of the Helmholtz equation was proposed, which combined the flexibility of the finite element method to deal with the boundary and region and the fast convergence of the spectral method. By means of the principle of weighted residuals, the basic theory and mathematical formulae of the Chebyshev spectral element method for transmission eigenvalue problems were obtained. The original problem was transformed into quadratic eigenvalue problems. Furthermore, several numerical examples were given to illustrate the effectiveness of the proposed method.

Key words: transmission eigenvalue problem; quadratic eigenvalue problem; spectral element method; Chebyshev basis function

作者简介：　戴海(1992—)，男，硕士生(通讯作者. E-mail: 1670112042@qq.com);潘文峰(1964—)，男，教授(E-mail: pan@mail.whut.edu).

戴海，潘文峰. 谱元法求解Helmholtz方程透射特征值问题[J]. 应用数学和力学, 2018, 39(7): 833-840.

DAI Hai, PAN Wenfeng. A spectral element method for transmission eigenvalue problems of the Helmholtz equation[J]. Applied Mathematics and Mechanics, 2018, 39(7): 833-840.

谱元法求解Helmholtz方程透射特征值问题*

引 言

1 二维透射特征值问题数值方法

1.1 数值模型

1.2 标准单元内函数Chebyshev谱逼近

1.3 单元刚度矩阵的生成

1.4 总体矩阵的生成

2 数 值 分 析

3 总结和发展