基于多项式插值的有限差分法求解Helmholtz方程透射特征值问题

Solving the Transmission Eigenvalue Problem of Helmholtz Equation by Finite Difference Method Based on Polynomial Interpolation

摘要: 有限差分公式在无网格方法求解微分方程数值解中起着重要作用。本文针对Helmholtz方程透射特征值问题，通过多项式插值来创建有限差分公式。本文运用一种简单实用的节点分布，既保证多元多项式插值的唯一可解性，又使矩阵为三角矩阵，以便构造的基本多项式化为Lagrange基多项式。最后给出了外透射特征值问题的数值算例。

Abstract: Finite difference formulas play an important role in the numerical solution of differential equations using meshless methods. In this paper, aiming at the transmission eigenvalue problem of Helmholtz equation, a finite difference formula is created by polynomial interpolation. Then, a simple and practical node distribution is used, which not only guarantees the unique solvability of multivariate polynomial interpolation, but also makes the matrix a triangular matrix so that the basic polynomials constructed can be transformed into Lagrange basis polynomials. Finally, a numerical example solving the eigenvalue problem of external transmission is given.

1. 引言

无网格方法(Mesh-less method)是相对于有网格方法的称谓 [1]，即在数值计算中不需要生成网格，而是按照一些任意分布的坐标点构造插值函数离散控制方程，方便地模拟各种复杂形状的流场。该方法大致可分成两类：一类是以Lagrange方法为基础的粒子法，如光滑粒子流体动力学法 [2]，和在其基础上发展的运动粒子半隐式法 [3] 等；另一类是以Euler方法为基础的无格子法，如无格子Euler/N-S [4] 算法和无单元Galerkin法 [5] 等。由于这种方法不需要划分网格，灵活多变，可以方便地利用坐标点计算模拟复杂形状流场计算，自从它提出以来就得到了快速发展。

虽然大多数问题都得到了解决，但在高雷诺数流动时提高数值计算精度较困难，而且具体的实施过程值得商榷和改进。黄小为和吴传生 [6] 就无网格方法发展的一种基于数值微分方法，进一步构造Lagrange多项式插值函数，和其对求解高维问题提出的行之有效的改进手段 [7]。如果采用全局近视，其优点是容易编程以及保持谱精确性，但相应地所得的线性方程组其矩阵条件数是相当大的，矩阵是病态的；如果由给出的局部节点构造插值多项式，则牺牲了谱精确性 [8]。在此数值微分方法中，对单个变量的情形，基于Lagrange插值是很容易实现的 [9] [10] [11]；在高维情形，由于多个变量插值的复杂性，不同求解区域的插值的唯一性 [12]，其唯一可解性是当且仅当线性方程组的范德蒙行列式不等零。这样要求插值的局部节点需要具有特殊的分布才能保持多变量插值的存在性和唯一性 [13] [14]。本研究主要是基于黄小为和吴传生求解Laplace方程以及热传导方程的工作 [7]，以及Gasca和Maeztu的工作 [15]，求解Helmholtz方程透射特征值问题。

散射问题和逆散射问题是数学物理研究领域中的经典问题，透射特征值问题源于非均匀的逆散射理论 [16]，是波场逆散射理论的重要内容。透射特征值问题最早是由Kirsch [17] 提出的，后来Colton和Monk [18] 对此做了更系统的研究。

本文结构如下：在第一节提出了一个简单且实用的节点分布并确定插值节点集；在第二节构造有限差分公式；在第三节给出外透射特征值问题的数值算例。

2. 节点分布和插值节点集

2.1. 节点分布

为了方便插值节点集的确定和插值多项式的构造，我们采用一种节点分配策略 [7]。在区域 $\Omega \subset R^2$ 上，为获得一组并行的线型节点，用 $\Delta$ 表示平面区域的离散化级别。为简单起见，我们假设 $\Omega =\left\{(x,y)|a_1\le x\le a_2,g_1\left(x\right)\le y\le g_2\left(x\right)\right\},a_1=\inf \left\{x|\left(x,y\right)\in \Omega \right\},a_2=\sup \left\{x|\left(x,y\right)\in \Omega \right\}$。

首先确定区域 $\Omega$ 在投影x轴上的节点数

$I=\{\begin{array}{l}1,a_1=a_2,\\ 2,0<a_2-a_1<\Delta ,\\ \left[\frac{a_2-a_1}{\Delta }\right]+1,a_2-a_1\ge \Delta ;\end{array}$

并在投影间隔上生成等距节点

$x_i=\{\begin{array}{l}{a}_1,I=1,\\ a_1+\frac{a_2-a_1}{I-1}\left(i-1\right),I>1,\end{array}i=1,\cdots ,I.$

类似的，我们生成一系列通过这些节点并平行于y轴的直线，并在直线和 $\Omega$ 上生成等距节点：

$J_i=\{\begin{array}{l}1,b_{i1}=b_{i2},\\ 2,0<b_{i2}-b_{i1}<\Delta ,\\ \left[\frac{b_{i2}-b_{i1}}{\Delta }\right]+1,b_{i2}-b_{i1}\ge \Delta ,\end{array}i=1,\cdots ,I;$

$y_{ij}=\{\begin{array}{l}{b}_{i1},J_i=1,\\ b_{i1}+\frac{b_{i2}-b_{i1}}{J_i-1}\left(j-1\right),J_i>1,\end{array}j=1,\cdots ,J_i.$

其中 $b_{i1}=\inf \left\{y|\left(x_i,y\right)\in \Omega \right\},b_{i2}=\sup \left\{y|\left(x_i,y\right)\in \Omega \right\}$。

最终，生成区域 $\Omega \subset R^2$ 上节点

$T={\displaystyle \underset{i=1}{\overset{I}{\cup }}\left\{\left(x_i,y_{ij}\right)|j=1,\cdots ,J_i\right\}.}$

在MATLAB中实现，得到如图1所示的圆形区域的节点分布

Figure 1. Node distribution

图1. 生成节点

共生成了1272个节点，执行时间约为0.013 s。

2.2. 插值节点集

在径向基方法中，通常选择最接近的一个邻居。类似的，在本小节中，对于给定的 $p\in \Omega$ 和 $n\in N$，我们着重于选择节点p附近的插值节点集 ${\chi }_n\left(p\right)$ [7]。

在1.1节中，生成了节点 $T=\left\{\left(x_i,y_{ij}\right)|i=1,\cdots ,I;j=1,\cdots ,J_i\right\}\subset \Omega$，然后将 ${\left\{J_i\right\}}_{i=1}^I$

降序排列为

$J_{m_1}\ge J_{m_2}\ge \cdots \ge J_{m_I}.$

为了保证 ${\chi }_n\left(p\right)$ 存在，让 $\Delta >0$ 为足够小的常数，使得 $I\ge n+1$ 并且

$J_{m_k}\ge n+2-k,k=1,\cdots ,n+1.$

为确定插值节点集，给定降序排列 $r=n+1,n,\cdots ,1$，先令集合 $S=\left\{1,\cdots ,I\right\}$，后构造 $V=\left\{i\in S|J_i\ge r\right\}$，然后通过下述规则迭代

$k\left(i\right)=\underset{1\le l\le J_i-r+1}{\arg \min }\left|v-\frac{y_{i,l}+y_{i,l+r-1}}{2}\right|,i\in V;$

$i_r=\underset{i\in V}{\arg \min }\left({\left(u-x_i\right)}^2+{\left(v-\frac{y_{i,k\left(i\right)}+y_{i,k\left(t\right)+r-1}}{2}\right)}^2\right);$

${\rho }_r\left(p\right)=\left\{\left(u_{i_r},v_{i_r,s}\right)|s=k\left(i_r\right),\cdots ,k\left(i_r\right)+r-1\right\};$

$S=S\backslash \left\{i_r\right\}.$

最终获得插值节点集

${\chi }_n\left(p\right)={\displaystyle \underset{r=1}{\overset{n+1}{\cup }}{\rho }_r\left(p\right)}=\left\{\left(u_r,v_{rs}\right)|\left(r,s\right)\in I_n^2\right\}\subseteq T.$

在MATLAB中实现，图2给出了圆域中给定点的插值节点，

Figure 2. Interpolation node set, denotes the given point p = (0.08, 0.45)

图2. 插值节点集p = (0.08, 0.45)

给定节点 $p=\left(0.08,0.45\right),n=6$，共生成了28个插值节点，执行时间约为0.011 s。

3. 有限差分公式

给定插值节点集

${\chi }_n\left(p\right)=\left\{q_{\beta }=\left(u_r,v_{rs}\right)|\beta =\left(r,s\right)\in I_n^2\right\},$

对于 $\alpha =\left(i,j\right)\in I_n^2$，令

${\xi }_i\left(x\right)=\{\begin{array}{l}1,i=1\\ {\displaystyle \underset{r=1}{\overset{i-1}{\prod }}\frac{x-u_r}{u_i-u_r},i>1,}\end{array}$ (1)

${\eta }_{ij}\left(y\right)=\{\begin{array}{l}1,j=1\\ {\displaystyle \underset{s=1}{\overset{j-1}{\prod }}\frac{y-v_{is}}{v_{ij}-v_{is}},j>1,}\end{array}$ (2)

关于 ${\chi }_n\left(p\right)$ 的插值空间 ${\Pi }_n^3$，可以用下面的牛顿基多项式来表示

${\phi }_{\alpha }\left(x,y\right)={\xi }_i\left(x\right){\eta }_{ij}\left(y\right),\alpha =\left(i,j\right)\in I_n^3.$ (3)

很容易证明

${\phi }_{\alpha }\left(q_{\beta }\right)=\{\begin{array}{l}1,\beta =\alpha \\ 0,\beta \prec \alpha \end{array}\text{,}\alpha ,\beta \in I_n^2.$ (4)

构造关于 ${\chi }_n\left(p\right)$ 的拉格朗日基多项式

$l_{\alpha }\left(q_{\beta }\right)=\{\begin{array}{l}1,\alpha =\beta ;\\ \text{0,}\alpha \ne \beta \text{.}\end{array}$ (5)

我们引入列向量

$\Phi ={\left[{\phi }_{\alpha }\right]}_{\alpha \in I_n^2},L={\left[l_{\alpha }\right]}_{\alpha \in I_n^2},$

和方阵

$A={\left[\Phi \left(q_{\beta }\right)\right]}_{\beta \in I_n^2}.$

由(4)，A是单位上三角矩阵且 $\det A=1$。很明显

$\Phi =AL,$

即

$L=A^{-1}\Phi .$ (6)

给定 $f:R^2\to R$，关于 ${\chi }_n\left(p\right)$，f的拉格朗日插值多项式为

${\displaystyle \underset{\alpha \in I_n^2}{\sum }f\left(q_{\alpha }\right)l_{\alpha }\left(q\right),q\in R^2}.$ (7)

对(7)中的多项式微分，并代入 $q=p$，可得出f在点处的近似导数为

$f^{\left(\gamma \right)}\left(p\right)\approx {\displaystyle \underset{\alpha \in I_n^2}{\sum }f\left(q_{\alpha }\right)D^{\gamma }{l}_{\alpha }\left(p\right),}$ (8)

其中 $\gamma \in N_0^2,\left|\gamma \right|\le n$， $D^{\gamma }{l}_{\alpha }\left(p\right)$ 是 $l_{\alpha }\left(q\right)$ 的微分。

4. 数值实验

考虑在 $D\in R^2$ 上的外透射特征值问题

$\{\begin{array}{l}\Delta w+k^{2}n\left(x\right)w=0,\text{in}D,\\ \Delta v+k^{2}v=0,\text{in}D,\\ w-v=f,\text{on}\partial D,\\ \frac{\partial w}{\partial \gamma }-\frac{\partial v}{\partial \gamma }=g,\text{on}\partial D.\end{array}$ (9)

其中 $\gamma$ 为边界上的单位法向量， $n\left(x\right)$ 是折射率，非零常数k为传输特征值。

将(8)代入(9)，可以得到

$\left(\begin{array}{c}{L}_w\\ L_v\\ F\\ B\end{array}\right)\left(\begin{array}{c}{w}_0\\ w_1\\ v_0\\ v_1\end{array}\right)+\left(\begin{array}{l}{k}^{2}n\left(x\right)\\ k^2\\ 0\\ \text{0}\end{array}\right)\left(\begin{array}{c}{w}_0\\ w_1\\ v_0\\ v_1\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ f\\ g\end{array}\right).$ (10)

其中 $w_0,w_1,v_0,v_1$ 为介质1、2内点和边界点的未知向量， $N_{w_0},N_{v_0}$ 为两个介质内点个数， $N_{w_1},N_{v_1}$ 为边界点个数， $N_{w_1}=N_{v_1}$， $N=N_{w_0}+N_{w_1}+N_{v_0}+N_{v_1}$。离散拉普拉斯算子 $L_{w_{N_{w_0}\times N}}=\left(L_{w_0},L_{w_1},0,0\right),L_{v_{N_{v_0}\times N}}=\left(0,0,L_{v_0},L_{v_1}\right)$ 和 $B=\left(B_{w_0},B_{w_1},B_{v_0},B_{v_1}\right)$ 的离散逼近分为四个块，分别对应于介质1的内点和边界点和介质2的内点和边界点。同样， $F=\left(0,E_{N_{w_1}},0,-E_{N_{v_1}}\right)$。

因此(10)也可写成下述矩阵表达

$\left(\begin{array}{cccc}{L}_{w_0}+k^{2}n\left(x\right)& L_{w_1}+k^{2}n\left(x\right)& 0& 0\\ 0& 0& L_{v_0}+k^2& L_{v_1}+k^2\\ 0& E_{N_{w_1}}& 0& -E_{N_{v_1}}\\ B_{w_0}& B_{w_1}& B_{v_0}& B_{v_1}\end{array}\right)\left(\begin{array}{c}{w}_0\\ w_1\\ v_0\\ v_1\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ f\\ g\end{array}\right),$

即

$\left(\begin{array}{c}{w}_0\\ w_1\\ v_0\\ v_1\end{array}\right)={\left(\begin{array}{cccc}{L}_{w_0}+k^{2}n\left(x\right)& L_{w_1}+k^{2}n\left(x\right)& 0& 0\\ 0& 0& L_{v_0}+k^2& L_{v_1}+k^2\\ 0& E_{N_{w_1}}& 0& -E_{N_{v_1}}\\ B_{w_0}& B_{w_1}& B_{v_0}& B_{v_1}\end{array}\right)}^{-1}\left(\begin{array}{c}0\\ 0\\ f\\ g\end{array}\right).$ (11)

给定 $n=6,\Delta =\frac{1}{20},D=\left\{\left(x,y\right)|x^2+y^2\le 4\right\},\partial D=\left\{\left(x,y\right)|x^2+y^2=1\right\}$， $f={\text{e}}^x,g=1$ 代入MATLAB，得图3。

Figure 3. Numerical solution of transmission eigenvalues (three-dimensional)

图3. 透射特征值数值解(三维)

共生成5159个节点，用时173.2355 s。

5. 误差分析

根据黄小为的假设 [7]：

存在一个数 $\lambda \ge \text{1}$，任意两点 $p_1,p_2\in \Omega$ 可通过长度为 $\left|\Gamma \right|\le \Vert p_1-p_2\Vert$ 的可校正曲线Γ连接。

$\Omega \subset R^2$ 是满足此假设的紧致域，T是通过上述算法生成的一个节点集，给定 $n\in N,p\in \Omega ,{\chi }_n\left(p\right)={\left\{q_{\alpha }\right\}}_{\alpha \in I_n^2}$ 是最终生成的插值节点集， ${\left\{l_{\alpha }\right\}}_{\alpha \in I_n^2}$ 是 ${\chi }_n\left(p\right)$ 生成的拉格朗日基。又 $\gamma =\left({\gamma }_1,{\gamma }_2\right)\in N_0^2,\left|\gamma \right|\le n$，若 $p\in \text{Rect}\left({\chi }_n\left(p\right)\right)\cap \Omega ,f\in C^{n,1}\left(\Omega \right)$，那么当

$C\left(\tau ,\kappa ,\gamma ,\lambda \right)=\frac{2\frac{3n+1}{2}{\lambda }^n}{\left(n-1\right)!}{\tau }^{n+1-\left|\gamma \right|}{C}_1\left(\tau ,\kappa ,\gamma \right)$ 时，

$\left|f^{\left(\gamma \right)}\left(p\right)-{\displaystyle \underset{\alpha \in I_n^2}{\sum }f\left(q_{\alpha }\right)l_{\alpha }^{\left(\gamma \right)}\left(p\right)}\right|\le C\left(\tau ,\kappa ,\gamma ,\lambda \right){\left|f\right|}_{n,1}{\Delta }^{n+1-\left|\gamma \right|},$

其中 $C_1\left(\tau ,\kappa ,\gamma \right)={\left(2\kappa \right)}^{\left|\gamma \right|}{\left({\tau }^n+1\right)}^{\left|I_n^2\right|-1}{\displaystyle \underset{\alpha =\left(i,j\right)\in I_n^2}{\sum }{\tau }^{\left|\alpha \right|-2}{T}_{i-1}^{\left({\gamma }_1\right)}\left(1\right)T_{j-1}^{\left({\gamma }_2\right)}\left(1\right)}$。

所以近似误差为 $\left|f^{\left(\gamma \right)}\left(p\right)-{\displaystyle \underset{\alpha \in I_n^2}{\sum }f\left(q_{\alpha }\right)D^{\gamma }{l}_{\alpha }\left(p\right)}\right|$。

6. 总结和发展

为了对PDE进行数值求解，本文提出了平行线型的节点分布，该方法便于选择插值节点集，并且保证了多元拉格朗日多项式插值的唯一可解性 [14] [15]。该插值节点集分布的最大优点就是，我们仅面对三角矩阵，就可通过构造的基础多项式获得拉格朗日基础多项式。算例表明，基于多项式插值的FD公式可以成功地应用于无网格法。

今后的工作可以从三维区域上的推广展开。