一种用于时间分数阶MIM方程的局部间断伽辽金方法：空间H1范数分析与参数识别

《Communications in Nonlinear Science and Numerical Simulation》：A Local Discontinuous Galerkin Method for the Time-Fractional MIM Equation: Spatial H1-Norm Analysis and Parameter Identification

　　**摘要**
在本文中，我们考虑了非线性时间分数阶移动/静止方程的数值解和参数识别问题，该方程常用于描述异质介质中的异常溶质传输。我们提出了一种高效的全离散数值方案，其中在空间上采用了局部不连续Galerkin方法，对于时间分数阶导数使用了分级网格上的平均L1近似，而对于一阶时间导数则采用了 первое порядка обратного Эйлера（first-order backward Euler）方法。对于初始时刻具有较弱规则性的解，我们证明了该方案在空间H1范数下的稳定性，并建立了最优误差估计。基于正向求解器的良好性能，我们设计了一种参数反演程序，该程序利用粒子群优化算法从观测数据中同时恢复分数阶数和传输系数。进行了大量的数值实验，包括受不同强度白噪声污染的合成示例以及经典的山区溪流示踪剂测试。数值结果表明，所提出的框架对于时间分数阶扩散过程的正向模拟和逆向建模都具有高精度和实用价值。

**引言**
理解和预测异质介质中的溶质传输对于地下水污染评估、碳氢化合物回收以及河流中的示踪剂分析至关重要。在这些系统中，观察到的突破曲线通常具有重尾的幂律特性，这无法用经典的平流-扩散方程再现。时间分数阶移动/静止（MIM）模型[1]源于具有幂律等待时间的连续时间随机游走，通过在被静止相中引入分数阶导数自然地捕捉了这种非菲克行为。因此，它已成为描述多孔含水层、河流系统和其他地球物理环境中异常传输的广泛使用的工具[2]。然而，由于初始时刻解的规则性较弱，MIM模型的准确数值模拟面临重大挑战，而从现场数据确定模型参数对于实际应用仍然是必要的。

在本文中，我们考虑了非线性时间分数阶MIM方程：
$$
\begin{cases}
\mathbf{u}_{t}(x, y, t) + \partial_{\lambda}\mathbf{u}_{t}(x, y, t) - \nabla\cdot(\mathbf{a}(x, y)\nabla\mathbf{u}_{t}(x, y, t) = \mathbf{f}(x, y, t), & (x, y, t) \in \Omega \times \mathbf{P}, \\
\mathbf{u}_{t}(x, y, 0) = \mathbf{u}_0(x, y), & (x, y) \in \Omega, \\
\mathbf{u}_{t}(x, y, t) = 0, & (x, y) \in \partial\Omega \times \overline{\Omega},
\end{cases}
$$
其中 $\alpha \in (0, 1)$，$\mathbf{P} = (0, \tau]$ 且 $\tau > 0$，$\Omega \subset \mathbb{R}^2$ 是一个有界的凸多边形域，$\overline{\Omega} = \Omega \cup \partial\Omega$。系数 $\mathbf{a}(x, y)$ 是光滑的，并满足 $0 < \mathbf{a}_\infty \leq \mathbf{a}(x, y) \leq \bar{\mathbf{a}$。Caputo分数阶导数 $\partial_{\lambda}\mathbf{u}_{t}$ 定义为
$$
\partial_{\lambda}\mathbf{u}_{t}(x, y, t) = \frac{1}{\Gamma}(1 - \lambda)\int_{t_0}^t \partial_{\sigma}\mathbf{u}_{t}(x, y, s)\partial_{\sigma}s(t - s) - \lambda\int_{s}^t \mathbf{u}_{s}, \quad 0 < \lambda < 1.
$$
非线性项 $\mathbf{f}$ 满足Lipschitz条件：$|\mathbf{f}(u_1) - \mathbf{f}(u_2)| \leq L|u_1 - u_2|$，其中 $L \geq 0$ 是Lipschitz常数。

时间分数阶MIM方程的一个主要数值挑战是初始时刻解的规则性较弱。基于Zheng等人的分析[3]，我们假设解具有以下弱规则性：存在一个常数 $\sigma > 0$，使得对于所有 $t \in \mathbf{P}$，有
$$
\| \mathbf{u}_t(\cdot, t) \|_{H^2(\Omega)} \leq C_1, \quad \| \mathbf{u}_{t}^2(\cdot, t) \|_{H^2(\Omega)} + \| \partial_{\sigma}\mathbf{u}_t(\cdot, t) \|_{H^2(\Omega)} + \| \mathbf{u}_{t}^2(\cdot, t) \|_{H^2(\Omega)} \leq C_2,
$$
这种弱规则性行为直接促使我们在分级时间网格上使用离散方案来近似时间分数阶导数，从而实现最优的收敛阶数。

分数阶MIM模型的数值模拟已经得到了广泛研究。例如，Qiao和Cheng[2]开发了一种快速的Crank-Nicolson L1有限差分方案，Jiang[4]构造了一种交替方向隐式紧凑差分方法，Yu和Chen[5]提出了一种用于分数阶移动/静止Allen-Cahn方程的平均L1公式，Salama和Fairag[6]采用了一种显式解耦群方法，Golbabai[7]引入了一种无网格径向基函数技术，而Fardi和Ghasemi[8]将时间有限差分与空间谱方法结合起来。另一方面，为了有效处理时间分数阶偏微分方程固有的初始奇异性，非均匀时间离散化已被证明是有效的方法。在这方面，有几项先前的工作。Stynes等人[9]对分级网格上的有限差分方案进行了严格的误差分析，为后续研究提供了重要基础，Liao等人[10]通过建立离散分数阶Gr?nwall不等式和进行全局一致性分析，进一步丰富了理论框架，从而获得了更准确的误差估计。在这些时间离散化策略的基础上，Jing和Li[11]提出了基于块中心有限差分和两种温度调整的L1离散化的数值方法，用于具有时间依赖系数的温度调整的次扩散方程，并证明了在非均匀网格上的稳定性和最优收敛率，而Li等人[12]分析了具有非光滑数据的高阶加权移位Grünwald-Letnikov方案，并通过引入校正项和积分-微分方法证明了修改后的方案在时间上实现了从二阶到四阶的最优收敛阶数。

局部不连续Galerkin（LDG）方法由Cockburn和Shu[13]提出，它基于将高阶方程重写为了一阶系统，然后使用标准DG框架进行离散化。这种方法自然地适应了悬挂节点、不规则元素几何形状和hp自适应性，同时保持了高阶精度和局部守恒性。LDG方案已成功应用于广泛的分数阶方程。Deng和Hesthaven[14]为分数扩散方程开发了一种惩罚性LDG方法，当分数阶数接近经典极限时，有效恢复了最优收敛阶数。对于非局部积分分数阶拉普拉斯方程，Nie和Deng[15]为二维问题提出了一个LDG框架，提供了完整的稳定性分析和分数阶的收敛估计。在非线性模型的背景下，Li和Liu[16]设计并分析了一种半离散LDG方案用于空间分数阶Allen-Cahn方程，为光滑解建立了严格的误差估计。值得注意的贡献还包括Li等人[17][18]，他们通过将空间LDG离散化与非均匀L1时间步进相结合，为一维和二维时间分数阶Burgers方程制定了稳定且收敛的完全离散框架。为了解决分数阶Rayleigh-Stokes问题中固有的初始奇异性，Agrawal等人[19]采用了一种最优收敛的LDG方法并确立了其适定性。将这种方法扩展到时间分数阶波传播，Li等人[20]为具有分数衰减的超声传播方程构建了两种新颖的半离散和完全离散的LDG方案，获得了最优收敛率。此外，Zhang等人[21]为二维非线性分数扩散方程开发了高阶LDG方法与二阶时间步进方案相结合，并提供了严格的先验误差估计。Li和Hou[22]还引入了在分级网格上的快速L1近似与LDG方法的结合，实现了时间分数阶双曲系统的最优收敛性和稳定性。

上述工作展示了LDG框架在分数阶模型中的良好适用性。然而，当将该框架应用于分级时间网格上的时间分数阶MIM方程时，直接与经典的L1或L2-1σ近似相结合会遇到结构困难。这些经典方案产生的离散卷积系数在非均匀网格上是正的但本质上是非对称的。在LDG框架中，全局能量估计需要一个对称的正定双线性形式来通过数值通量吸收时空交叉项。经典系数的不对称性阻止了这种闭合。因此，许多现有的LDG分析要么将其误差估计限制在L2范数内，要么仅在均匀时间网格上有效。在分级网格上，它们通常依赖于人为强的时间规则性假设来补偿这种结构不匹配[23][24]。

本文的主要创新在于扩展了LDG框架的适用性。通过采用平均L1公式[25][26][27]，恢复了所需的对称性。主要的理论贡献是为变量系数MIM方程开发了一个专门的数值分析框架。通过结合高斯-拉道投影算子和特定的测试函数，成功地将平均L1核的对称性纳入了L2能量论证中。此外，证明中的一个关键步骤是将光滑系数 $\mathbf{a}(x, y)$ 分解为常数部分和几何上小的扰动。这种分解，加上对称的核，使我们能够控制由LDG数值通量引入的边界跳跃项。因此，在初始时刻的弱规则性条件（1.3）下，建立了无条件稳定性和最优的空间H1范数误差估计，而无需引入任何人为光滑性假设。

除了对正向问题的数值分析外，该方案的实用性还通过将其与参数识别框架相结合得到了进一步的证明。优化算法（如Levenberg-Marquardt（LM）方法和粒子群优化（PSO）[28][29][30][31]的成功，严重依赖于正向求解器提供的梯度质量。基于均匀网格的传统求解器通常无法准确解决初始的奇异性，导致梯度信息失真，从而削弱了反演程序的可靠性。相比之下，所提出的平均L1-LDG求解器即使在规则性较弱的条件下也能提供高精度的场近似，从而为优化算法提供了稳定和可靠的梯度。通过将这种高精度正向求解器与LM和PSO算法相结合，构建了一个强大的反演框架，能够成功估计分数阶数和扩散系数。通过涉及噪声污染的合成数据和真实示踪剂测量的数值实验，验证了该反演框架的准确性和可靠性，证实了其在实际应用中表征异常扩散的潜在能力。

本文的其余部分安排如下。第2节介绍预备知识。第3节描述了分级网格上的平均L1-LDG方案，分析了其在空间H1范数下的稳定性和误差估计，并提供了数值示例来验证理论结果。第4节将平均L1-LDG方法应用于逆问题。使用LM和PSO算法从观测数据中识别分数阶数和传输系数，进一步证明了该数值方案的工程实用性。第5节给出了结论。

**预备知识**
我们首先定义了后续分析所需的投影和插值算子。设计算域为矩形 $\overline{\Omega} = [a, b] \times [c, d]$，它被划分为由 $N = N_x \times N_y$ 个矩形元素 $\mathbf{E}_{ij} = [x_{ij-1/2}, x_{ij+1/2} \times [y_{j-1/2}, y_{j+1/2}] = \mathbf{I}_{ij}$, 对于 $i = 1, \ldots, N_x$ 和 $j = 1, \ldots, N_y$。网格节点满足 $a = x_{1/2} < x_{3/2} < \ldots < x_{N/2} = b$ 和 $c = y_{1/2} < y_{3/2} < \ldots < y_{N/2} = d$。局部网格宽度为 $\Delta_x_{ij} = x_{ij+1/2} - x_{ij-1/2}$ 和 $\Delta_y_{j} = y_{j+1/2} - y_{j-1/2}$，$\hat{N} = \max\{ \Delta_x_{ij}, \Delta_y_{j} \}$ 表示最大边长。

**数值方案**
在介绍我们改进的方法之前，我们简要回顾了分级网格上的标准L1-LDG离散化，这将作为基线，以突出时间精度的提升。

**所提出数值方法（3.19）的应用**
在许多实际和工程设置中，控制方程中的某些参数或组分（如源项、边界数据或初始条件）通常是事先不可知的，必须被视为未知的。因此，开发有效的计算策略来识别未知参数或从实验或观测数据中重建未知功能形式引起了相当大的关注，并构成了一个重要的研究课题。

**结论**
在本文中，我们构建并严格分析了分级时间网格上的平均L1-LDG方法，用于非线性时间分数阶移动/静止（MIM）方程。通过利用平均L1离散算子的对称正定性，我们在没有人为规则性条件的情况下建立了H1范数下的最优误差估计。我们克服的关键困难是经典的L1和L2-1σ公式仅产生

**作者贡献声明**
王玉杰：写作——原始草稿、软件、方法论。王震：写作——审稿与编辑、方法论、资金获取。陈俊玉：软件。

**利益冲突声明**
作者声明他们没有已知的竞争财务利益或个人关系可能影响本文报道的工作。

**学校信息**
王玉杰 | 王震 | 陈俊玉
江苏大学数学科学学院，镇江，212013，中国