基于稀疏网格的高维有限元载荷向量数值积分：带振荡核函数的自适应节点加密策略
题目描述
考虑计算一个出现在三维波动方程有限元离散化中的高维积分。该积分具有振荡核函数，其形式为：
\[I = \int_{\Omega} f(\mathbf{x}) \, e^{i \kappa \|\mathbf{x} - \mathbf{x}_0\|} \, d\mathbf{x} \]
其中：
- \(\Omega \subset \mathbb{R}^3\) 是一个三维有限元单元（例如，一个六面体或四面体）。
- \(f(\mathbf{x})\) 是一个光滑但可能具有边界层特征的标量函数，例如一个多项式形状函数。
- \(e^{i \kappa \|\mathbf{x} - \mathbf{x}_0\|}\) 是振荡核，\(\kappa\) 是一个较大的波数（即高频振荡），\(\mathbf{x}_0\) 是空间中的一个固定点。
- 该积分在计算电磁散射或声学问题的有限元矩阵元素（即载荷向量项）时出现。
我们的目标是高效、准确地计算该积分。难点在于：
- 高维性：三维积分，传统张量积高斯求积节点数随维度指数增长（维度灾难）。
- 高频振荡：当 \(\kappa\) 很大时，被积函数剧烈振荡，标准高斯求积需要大量节点才能捕捉振荡，效率极低。
- 局部特征：\(f(\mathbf{x})\) 可能在单元边界附近变化剧烈（边界层），需要局部细化。
要求设计一个基于稀疏网格的自适应数值积分策略，结合振荡核的特性，通过自适应加密节点来平衡精度和计算成本。
解题过程
步骤1: 问题转化与基准方法分析
首先，将问题映射到标准参考单元。对于六面体单元，通常通过等参变换 \(\mathbf{x} = \mathbf{G}(\boldsymbol{\xi})\) 将物理单元 \(\Omega\) 映射到参考立方体 \(\hat{\Omega} = [-1, 1]^3\)。积分变为：
\[I = \int_{\hat{\Omega}} f(\mathbf{G}(\boldsymbol{\xi})) \, e^{i \kappa \|\mathbf{G}(\boldsymbol{\xi}) - \mathbf{x}_0\|} \, |J(\boldsymbol{\xi})| \, d\boldsymbol{\xi} \]
其中 \(J\) 是变换的雅可比行列式。记 \(F(\boldsymbol{\xi}) = f(\mathbf{G}(\boldsymbol{\xi})) |J(\boldsymbol{\xi})|\)，\(\phi(\boldsymbol{\xi}) = \kappa \|\mathbf{G}(\boldsymbol{\xi}) - \mathbf{x}_0\|\)，则积分简化为：
\[I = \int_{[-1,1]^3} F(\boldsymbol{\xi}) \, e^{i \phi(\boldsymbol{\xi})} \, d\boldsymbol{\xi} \]
这是一个标准的三维高振荡积分。
基准方法对比：
- 张量积高斯求积：在每一维使用 \(n\) 点高斯-勒让德规则，总节点数为 \(n^3\)。当 \(n\) 需要很大以解析振荡时，计算量爆炸。
- 标准稀疏网格求积（Smolyak算法）：利用一维求积规则的张量积的稀疏组合，节点数从 \(O(n^3)\) 降为 \(O(n (\log n)^{2})\)。但对于高频振荡函数，均匀的稀疏网格可能仍然需要很高的精度等级 \(n\)，导致节点数过多。
因此，我们需要一个能针对振荡函数和光滑函数 \(F\) 的边界层进行自适应局部加密的稀疏网格方法。
步骤2: 自适应稀疏网格求积框架
我们采用基于嵌套规则的稀疏网格构造，例如 Clenshaw-Curtis 规则或高斯-帕特森规则，它们的节点集是嵌套的（即低精度等级的节点集是高精度等级节点集的子集）。这便于自适应细化。
设有一维嵌套求积规则，其精度等级为 \(l\)（\(l=1,2,\dots\)），对应的节点集为 \(X_l\)，权重为 \(w_l\)。规则越精确，\(l\) 越大。
对于 \(d=3\) 维，稀疏网格求积公式可以写为：
\[Q_{L}^d f = \sum_{l_1+\cdots+l_d \leq L} \Delta_{l_1} \otimes \cdots \otimes \Delta_{l_d} f \]
其中 \(L\) 是总精度水平，\(\Delta_l = Q_l - Q_{l-1}\) 是第 \(l\) 层与第 \(l-1\) 层规则的差值算子（\(Q_0 = 0\)）。这个公式是 Smolyak 算法的另一种表述。它只组合那些各维精度等级之和不超过 \(L\) 的张量积差值，从而大幅减少节点。
关键思想：我们不预先固定总水平 \(L\)，而是从一个小网格开始，然后根据局部误差指示器自适应地添加那些对减少误差贡献最大的子网格（即 \(\Delta_{l_1} \otimes \cdots \otimes \Delta_{l_d}\) 对应的节点集）。
步骤3: 针对振荡核的误差指示器设计
标准的自适应稀疏网格通常使用函数值的差值或高阶导数作为误差指示。但对于高频振荡函数 \(e^{i \phi(\boldsymbol{\xi})}\)，其导数很大，直接使用会导致过度细化。
我们需要一个更能反映积分误差的指示器。考虑一个子网格对应的差值项 \(\Delta_{l_1} \otimes \Delta_{l_2} \otimes \Delta_{l_3} (F e^{i \phi})\)。其贡献的绝对值可以作为该子网格重要性的度量。但直接计算这个差值需要在该子网格的所有节点上计算被积函数，成本高。
我们利用振荡积分的特性进行简化：
- 相位函数近似：在局部区域，如果相位函数 \(\phi(\boldsymbol{\xi})\) 变化剧烈，其振荡主导了积分行为。我们可以用线性近似：\(\phi(\boldsymbol{\xi}) \approx \phi(\boldsymbol{\xi}_c) + \nabla \phi(\boldsymbol{\xi}_c) \cdot (\boldsymbol{\xi} - \boldsymbol{\xi}_c)\)，其中 \(\boldsymbol{\xi}_c\) 是该子网格对应区域的中心。
- 分离振荡部分：在该子网格区域，将积分视为 \(\int F(\boldsymbol{\xi}) e^{i \phi(\boldsymbol{\xi})} d\boldsymbol{\xi} \approx e^{i \phi(\boldsymbol{\xi}_c)} \int F(\boldsymbol{\xi}) e^{i \nabla \phi(\boldsymbol{\xi}_c) \cdot (\boldsymbol{\xi} - \boldsymbol{\xi}_c)} d\boldsymbol{\xi}\)。后半部分是一个带线性相位的积分。
- 误差指示器构造：对于给定的子网格（由多维指标 \(\mathbf{l} = (l_1, l_2, l_3)\) 标识），我们定义其误差指示器 \(\eta_{\mathbf{l}}\) 为：
\[ \eta_{\mathbf{l}} = \left| \Delta_{l_1} \otimes \Delta_{l_2} \otimes \Delta_{l_3} \left( F(\boldsymbol{\xi}) e^{i \nabla \phi(\boldsymbol{\xi}_c) \cdot \boldsymbol{\xi}} \right) \right| \]
注意这里我们移除了全局相位 $e^{i \phi(\boldsymbol{\xi}_c)}$（因为它不影响模长），并用线性相位 $\nabla \phi(\boldsymbol{\xi}_c) \cdot \boldsymbol{\xi}$ 代替了原相位。这样，差值运算作用于函数 $G(\boldsymbol{\xi}) = F(\boldsymbol{\xi}) e^{i \nabla \phi(\boldsymbol{\xi}_c) \cdot \boldsymbol{\xi}}$，该函数的振荡频率由局部波向量 $\nabla \phi(\boldsymbol{\xi}_c)$ 决定，在局部区域内近似为常数频率振荡。
这个指示器 \(\eta_{\mathbf{l}}\) 可以相对廉价地计算，因为它只涉及在子网格节点上计算 \(G(\boldsymbol{\xi})\)，而不需要精确的 \(e^{i \phi(\boldsymbol{\xi})}\)。
步骤4: 自适应算法流程
- 初始化：从最小的稀疏网格开始，例如 \(L=3\)（各维等级和不超过3）。计算当前网格的积分近似值 \(I_{\text{current}}\)。
- 计算候选子网格：找出所有当前网格的“子节点”子网格。即那些满足：其多维指标 \(\mathbf{l}\) 的每个分量都不小于当前网格中已有指标，且至少有一个分量比当前网格中对应指标大1，同时其他分量至少有一个与当前最大指标相同（这是稀疏网格的邻接细化关系）。为每个候选子网格计算误差指示器 \(\eta_{\mathbf{l}}\)。
- 选择与细化：将所有候选子网格按 \(\eta_{\mathbf{l}}\) 从大到小排序。选择前 \(k\) 个（例如 \(k=5\)）指示器最大的子网格，将它们加入到当前稀疏网格中。这意味着添加这些子网格对应的差值项 \(\Delta_{l_1} \otimes \Delta_{l_2} \otimes \Delta_{l_3}\) 到求积公式中。
- 更新积分值：更新积分近似：\(I_{\text{new}} = I_{\text{current}} + \sum_{\text{selected } \mathbf{l}} \Delta_{l_1} \otimes \Delta_{l_2} \otimes \Delta_{l_3} (F e^{i \phi})\)。注意，这里需要使用完整的被积函数 \(F e^{i \phi}\) 来计算这些新增的差值项。
- 误差估计与终止：可以估算全局误差为当前所有候选子网格的 \(\eta_{\mathbf{l}}\) 之和的某个比例，或者直接监测积分值的变化量 \(\| I_{\text{new}} - I_{\text{current}} \|\)。当该变化量小于预设容差 \(\epsilon\) 时，算法终止。否则，令 \(I_{\text{current}} = I_{\text{new}}\)，返回步骤2。
步骤5: 处理边界层（函数 \(F\) 的奇异性）
如果 \(F(\boldsymbol{\xi})\) 在单元边界（如 \(\xi_i = \pm 1\)）附近存在边界层（急剧变化），标准稀疏网格在边界处的分辨率可能不足。我们的策略可以自然地处理这一点：
- 在边界层区域，函数 \(F\) 变化剧烈，其导数很大。即使没有振荡核，误差指示器 \(\eta_{\mathbf{l}}\) 中由 \(F\) 贡献的部分也会变大（因为 \(G(\boldsymbol{\xi})\) 的振幅变化大）。
- 因此，算法会自动地在边界层方向（即对应于边界法向的维度）上选择更高等级的子网格进行加密，从而在边界附近布置更多节点。
步骤6: 总结与优势
本方法的核心是将振荡核的局部相位线性化与稀疏网格的自适应细化相结合。
- 优势1：通过稀疏网格，避免了维度灾难，节点数远少于张量积网格。
- 优势2：通过基于局部线性相位近似的误差指示器，能够识别出那些对高频振荡积分贡献最大的区域（通常是相位变化最快的方向），并进行针对性加密。
- 优势3：自适应机制同时照顾了光滑函数 \(F\) 的边界层特征，实现综合的局部节点优化布局。
- 结果：对于一个给定的精度要求，本方法能够使用比均匀稀疏网格或张量积网格少得多的函数求值次数，高效计算出带振荡核的高维有限元载荷向量积分。
通过这种方式，我们为高频波动问题有限元分析中的核心积分计算提供了一个实用且高效的数值工具。