基于FFT加速的基本解法计算非均匀外部亥姆霍兹问题的狄利克雷-诺伊曼映射

《Mathematics and Computers in Simulation》：FFT-accelerated computation of the Dirichlet-to-Neumann map for inhomogeneous exterior Helmholtz problems using the method of fundamental solutions

　　

摘要

本文提出了一种高效的数值方法，用于计算具有非齐次项的二维Helmholtz方程外Dirichlet问题所对应的Dirichlet到Neumann（DtN）映射。通过基本解方法（MFS）对方位人工边界外部子域内的解进行逼近。当源点和配点均匀分布在同心圆上时，离散化产生的系数矩阵成为循环矩阵，这使得可以通过快速傅里叶变换（FFT）高效地计算离散DtN映射。通过将DtN映射定义的透明边界条件（TBC）引入有限元格式，原始定义在无界域上的外问题被归结为定义在有界计算域上的等价问题，进而可以用有限元法（FEM）求解。数值算例表明：所提出的基于MFS的方法能够高精度地逼近DtN映射；基于FFT的算法能够高效地计算所得的离散DtN映射；将DtN映射定义的TBC纳入有限元方法后，能够获得无界域上具有非齐次项的外Helmholtz问题的精确解。

引言

椭圆型偏微分方程的外边界值问题在诸多应用中有所涉及，包括波动传播和散射。一个典型例子是外非齐次Helmholtz问题，该问题出现在当局部源位于障碍物附近的时谐声波辐射中。例如，由振动设备或扬声器在障碍物附近发出的波可以由Helmholtz方程描述，其外部域中包含一个紧支集的非齐次项。由于辐射场传播到无界域中，对辐射条件的准确处理至关重要。

在偏微分方程的数值方法中，有限元法（FEM）是最可靠且最广泛使用的方法之一。然而，FEM天然地定义在有界计算域上，因此不能直接应用于无界域上的外问题。因此，为了用FEM求解外Helmholtz问题，有必要将原始问题替换为定义在有界域上的等价问题。

在此背景下，Dirichlet到Neumann（DtN）映射起着核心作用。DtN映射将施加在人工边界上的Dirichlet数据与该边界上相应的法向导数关联起来。一旦在人工边界上获得DtN映射，原始外问题就可以被重新表述为定义在有界域上、并在人工边界上具有透明边界条件（TBC）的等价边界值问题 [1], [2]。这使得可以在有界计算域上使用FEM来求解无界域上的原始问题。

外边界值问题的一种经典方法是Dirichlet-Neumann交替法（D-N方法）[3], [4], [5], [6]。这是一种迭代域分解方法 [7]，其中人工边界将计算域分为两个子域：非齐次方程在内部子域中求解，齐次方程在外部子域中求解。其收敛速度强烈依赖于用于更新人工边界上边界值的松弛参数。

对于二维Poisson方程，当原始问题的边界和人工边界都是同心圆时，产生最快收敛的最优常数松弛参数是已知的 [8]。作者进一步推导了一个最优可变松弛参数，其中松弛参数在每次迭代中被最优地选取 [5]。虽然圆形边界的假设简化了理论分析，但一般而言原始问题的边界并不一定是圆形的。因此，考虑更一般形状的人工边界是有意义的 [2]。另一种可能的方法是使用保形映射将问题转化为具有圆形边界的问题 [9]，从而可以应用最优可变松弛参数。

然而，尽管保形映射对二维Poisson方程有效，但其对Helmholtz方程的应用变得困难，因为波数会被转化为空间变化的系数。此外，保形映射在三维中不可用。这些困难促使开发一种不依赖于局限于圆形边界的迭代松弛方案或保形映射的直接数值方法。

在本文中，我们开发了这样一种直接数值方法，用于求解二维非齐次Helmholtz方程的外Dirichlet问题。如上所述，通过引入一个圆形人工边界，我们通过DtN映射表示人工边界外部无界子域的影响，并将所得的条件作为TBC施加在人工边界上。这将原始外问题归结为有界域问题，进而可以用FEM求解。

为了数值构建DtN映射，我们采用基本解方法（MFS）[10], [11], [12], [13], [14], [15]，这是一种偏微分方程的无网格边界型方法。在MFS中，解被近似为基本解的有限线性组合，其奇点（称为源点）放置在计算域外部。所得的近似在整个域内精确满足控制方程。因此，展开系数仅通过在边界上的配点强制边界条件来确定。由于MFS不需要网格生成或体积积分，其实现简单，并且众所周知对于具有光滑几何和光滑解的问题能够提供高精度的近似。

在构建DtN映射中所采用的MFS属于更广泛的径向基函数（RBF）方法类。虽然RBF方法可以提供高精度且高效的离散化 [16]，并且在适当的节点配置下可能允许利用循环矩阵结构进行计算 [9], [17]，如下文所述，但本研究将基于MFS的DtN映射与FEM耦合，因为FEM为非齐次边界值问题和复杂几何提供了成熟的变分框架。这种耦合也表明，本文提出的计算DtN映射的方法可以自然地纳入标准计算框架中，尽管它本质上并不局限于FEM。完全的无网格公式将需要关于基函数、节点分布、形状参数和稳定化的额外选择。对这些问题的系统研究将构成一项独立的研究，超出了本文的范围。

对于圆形人工边界，精确的DtN映射也可以替代地用Fourier级数表示 [3]。这种解析构造非常有效，但它依赖于变量分离和外解的显式模态表示。边界积分方程也为构建DtN映射提供了标准框架 [4]，尽管其数值实现通常涉及奇异或近奇异的边界积分。

此处采用MFS是因为，如上所述，它在无需对外域进行网格划分或奇异积分的情况下，满足外Helmholtz方程和辐射条件。此外，它直接构建离散DtN映射，并在旋转对称性下允许基于FFT的加速，如下文所述。虽然它可能因源点位置的不同而遭受病态问题，但对本目的而言仍然适用。

MFS与快速傅里叶变换（FFT）的结合长期以来被用于在径向对称域中高效求解各种边界值问题。这类方法通常被称为矩阵分解算法（MDAs）；例如，见文献[18]第9章。在本研究中，采用MDA通过利用所得矩阵的循环结构来高效计算离散DtN映射 [9], [17], [19], [20]。

在作者之前的论文中，已在文献[21]中提出了齐次外Helmholtz问题的基于FFT的DtN映射数值计算的相关方法。在本文中，我们将这一研究方向扩展到非齐次Helmholtz方程。主要要点有两方面：第一，对MFS应用MDA来准确高效地在人工边界上构建DtN映射；第二，将所得的由DtN映射定义的TBC纳入有限元格式，从而使得无界域上的非齐次外Helmholtz问题可以在有界计算域中求解。

本文其余部分的组织如下。第2节建立二维非齐次Helmholtz方程的外边界值问题，并引入圆形人工边界和DtN映射。第3节给出MFS格式，用于近似人工边界外部的解，以及由此产生的DtN映射离散化。第4节开发基于FFT的算法，通过利用矩阵的循环结构来高效计算离散DtN映射。第5节描述具有由DtN映射定义的TBC的有限元离散化。第6节给出数值算例，展示了基于MFS的DtN映射的高精度、基于FFT计算的高效性、将DtN映射定义的TBC成功纳入齐次外问题的有限元格式，以及所得方法对非齐次外问题的适用性，包括一个具有物理动机的算例。最后，第7节总结全文。

节选片段

**问题建立**

在全文中，我们通过映射 $(x,y) \mapsto x+iy$ 将 $\mathbb{R}^2$ 与 $\mathbb{C}$ 等同，并记 $z = x+iy$，其中 $i := \sqrt{-1}$ 表示虚数单位。设 $R \in C[0,2\pi)$ 为给定函数。我们考虑简单闭曲线 $\Gamma := \{R(\theta)e^{i\theta} : 0 \leq \theta < 2\pi\}$ 及由 $\Gamma$ 包围的单连通有界域 $\Omega_0 \subset \mathbb{C}$。不失一般性，我们假设 $0 \in \Omega_0$。我们考虑以下非齐次Helmholtz方程的外Dirichlet问题，其中待定函数 $u \in H^1_{loc}(\mathbb{C} \setminus \overline{\Omega_0})$：

$$-\left(\Delta + \kappa^2\right)u = f \quad \text{在} \; \mathbb{C} \setminus \overline{\Omega_0} \; \text{中}$$

$$u = g \quad \text{在} \; \Gamma \; \text{上}$$

**解的离散化**

对于 $\zeta \in \mathbb{C}$，Helmholtz算子 $-(\Delta + \kappa^2)$ 的基本解 $\Phi_\kappa(\cdot, \zeta)$ 定义为

$$-\left(\Delta + \kappa^2\right)\Phi_\kappa(z, \zeta) = \delta(z-\zeta)$$

在分布的意义下成立，其给出

$$\Phi_\kappa(z, \zeta) := \frac{i}{4}H_0^{(1)}(\kappa|z-\zeta|), \quad z, \zeta \in \mathbb{C}, \quad z \neq \zeta$$

其中 $\delta(z-\zeta)$ 表示集中在点 $\zeta$ 处的Dirac delta分布，$H_\nu^{(1)}$ 表示 $\nu$ 阶第一类Hankel函数。

我们现在简要证明基本解 $\Phi_\kappa(\cdot, \zeta)$ 满足Sommerfeld辐射条件 (2.4c)。固定源点 $\zeta$，并写 $z = r\hat{z}$，其中 $r = |z|$ 且 $|\hat{z}|$

**使用FFT的高效计算**

在本节中，我们描述一种计算离散DtN映射的高效方法。关键观察是：当源点和配点都在圆上均匀分布时，来自MFS离散化的矩阵成为循环矩阵。这种结构使得可以通过FFT高效地计算离散DtN映射 [9], [17], [19], [20]。

虽然循环矩阵的性质是标准的，但我们在此陈述相关细节以固定全文使用的约定

**具有透明边界条件的有限元格式**

在本节中，我们简要说明如何将透明边界条件（TBC）纳入定义在 $\Omega$ 中的内部问题的有限元法（FEM）中。将Helmholtz方程 (2.6a) 乘以试验函数 $w$ 的复共轭 $\bar{w}$ 并应用Green公式，我们得到

$$\int_\Omega \nabla u \cdot \nabla \bar{w} \, d\Omega - \kappa^2 \int_\Omega u\bar{w} \, d\Omega - \int_{\Gamma_0} (\Lambda u)\bar{w} \, d\Gamma = \int_\Omega f\bar{w} \, d\Omega$$

对所有满足 $w|_\Gamma = 0$ 的 $w \in H^1(\Omega)$ 成立。设 $V_h \subset H^1(\Omega)$ 为由分段线性函数组成的有限元空间。设 $f_h$ 和 $g_h$ 为……

**离散DtN映射的精度**

对于给定的Dirichlet数据 $\lambda$，我们求解问题 (2.4)。我们首先考察MFS近似 $v^{(N)}$ 对外部子域 $\Omega_e$ 中精确解 $v$ 的精度，然后评估人工边界 $\Gamma_0$ 上相应的Neumann数据 $\mu$ 的精度，从而验证基于MFS的离散DtN映射的精度 [21]。圆形人工边界 $\Gamma_0$ 的半径取为 $R_0 = 1$，波数设为 $\kappa = 30$。精确解由以下表达式给出：

$$v(re^{i\theta}) := \frac{H_1^{(1)}(\kappa r)}{H_1^{(1)}(\kappa)}\cos\theta$$

**结论**

在本文中，我们提出了一种直接数值方法，基于域分解方法结合数值构建的DtN映射，用于求解无界域上二维非齐次Helmholtz方程的外Dirichlet问题。DtN映射通过MFS进行离散化，从而得到循环矩阵表示，这使得可以通过FFT高效地应用该映射。

通过将 DtN 映射定义的 TBC 纳入有限元离散化中，致谢：作者衷心感谢审稿人仔细阅读并提出了富有建设性和帮助性的意见，这些意见极大地提高了稿件的质量。重田貴美