一种用于微非均匀多孔介质中波动方程的高效Hermite谱求解器

《Communications in Nonlinear Science and Numerical Simulation》：An Efficient Hermite Spectral Solver for Wave Equations in Micro-Heterogeneous Porous Media

　　侯泽鑫 | 李灿 | 郭世民
西安工业大学应用数学系，中国陕西省西安市710054

摘要
我们提出了一种高效的海mite谱方法，用于求解描述微异质多孔介质中波传播的积分-微分偏微分方程的初值问题。利用温和分数阶微分的性质，首先将具有特定记忆核的原始积分-微分方程重新表述为等效的时分数阶偏微分方程。利用Faedo-Galerkin技术，我们证明了在无界域中重新表述的模型的适定性。对于数值离散化，空间导数采用海mite谱方法处理，结果表明得到的半离散方案是无条件稳定的且具有谱精确性的。在时间上，二阶项通过中心差分方案离散化，而温和的Caputo分数阶导数则通过加权移位Grünwald-Letnikov公式近似。数值实验不仅验证了理论稳定性和收敛性分析，还确认了完整离散化的二阶时间精度。此外，模拟显示了该模型展示的多种动态行为。

引言
多孔介质中的波传播在地球物理学、声学和材料科学中具有基础重要性[1]。在充满粘性流体的微异质多孔介质中，扩散关系受到孔隙流体流动中粘性边界层存在的强烈影响[2]。扩散弛豫产生了类似的扩散关系[3]。根据Gurevich-Lopatnikov理论，薄层随机介质中快波的扩散行为源于主要的衰减机制：在材料界面将快纵波转换为扩散慢波[4]。实证发现和严格的数学分析共同表明，分数幂扩散律的出现是时间域中半阶分数阶导数的表现。对于沿正x方向传播的波u(x, t)，对于高波频率ω，分数幂扩散律给出[3]，[5]
κ=2π/a·c∞·ω^(1/2)·(1?iω)/(1/τ^R)，
其中κ=2π/a表示复波数，c∞是声速的高频极限，a是一个取决于声速的高频和低频极限以及τR的常数。关系式(1)提供了方程
1/c∞·?^2u(x, t)/?t + ?u(x, t) ? a·q^(1/2)·?^2u(x, t) = 0，
q^(1/2)(t) = t^(?1/2)/π，
其中*表示时间上的拉普拉斯卷积算子，即(f·g)(t) := ∫_0^t f(t?s)·g(s)·ds。在具有多尺度记忆效应的介质中，一般的扩散律可以描述为[6]
κ=1/c∞·ω + a·ω^(β?iω)·α，α∈(0,1)，β>0，
这产生了向前的波动方程
1/c∞·?^2u(x, t) + ?u(x, t) ? a·q^α·β(t)·u(x, t) = 0，
其核函数为
q^α·β(t) = t^α?1·Γ(α)·e^(-βt)，α∈(0,1)，β>0，
这提供了积分-微分方程[7]
1/c∞·?^{2(?^t??^α·β(t)*)·(?^t+?^α·β(t)*)·u(x, t) = 0。
借助卷积的性质，方程(3)可以重写为[8]
1/c∞·?^{2}u_t(x, t) ? Δu(x, t) + (2a·q^α·β(t) + a^2·q^α·β(t)·q^α·β(t))·?^{2}u_t(x, t) = 0。
在本文中，我们不失一般性地考虑以下问题[6]
?^{2}u_t(x, t) ? Δu(x, t) + K(t)·?^{2}u_t(x, t) = f(x, t)，
x∈R^d，t∈(0, T)，
其中x=(x1,?, x_d)（d=1,2,3）是无界域上的多维变量，K(t)表示卷积核。在本文中，我们考虑K(t)=2a·q^α·β(t) + a^2·q^α·β(t)*q^α·β(t)，或K(t)=2a·q^α·β(t)，采用核函数(2)。
由于在不同的领域和应用背景下有不同的需求，已经建立并数值求解了不同的声波模型，包括地球科学[9]、[10]、生物工程[11]、岩石和土壤力学以及金融建模[12]。与经典声波模型相比，分数阶微分方程能够更精确地捕捉介质中的记忆效应和遗传性质。特别是在波传播问题中，分数幂律衰减波动方程成功地揭示了损耗介质中遵循频率依赖的幂律关系的衰减机制[13]。鉴于其建模优势，许多学者对分数波方程进行了广泛研究并取得了显著进展[13]、[14]、[15]、[16]、[17]、[18]。
虽然在某些特定假设下可以获得方程(4)的数值解和解析解，但这仍然具有很大挑战性。幸运的是，回想一下温和的Caputo分数阶导数的变体[19]
0_CDt^α·λ·u(t) = 1/Γ(1?α)·∫_0^t e^(-λ(t?s)·(t?s)·(-α)·u'(s)·ds，α∈(0,1)，λ>0，
可以直接计算得到(2a·q^α·β(t) + a^2·q^α·β(t))·?^{2}u_t(x, t) = 2a_0_CDt^(1?α,β)·u(t) + a^2_0_CDt^(1?2α,β)·u(t)。在下文中，我们考虑在无界域上的模型(4)，其核函数形式为K(t)=2a·q^α·β(t) + a^2(q^α·β·q^α·β)(t)。即我们研究
{?^{2}u_t(x, t) ? Δu(x, t) + 2a_0_CDt^(1?α,β)·u(t) + a^2_0_CDt^(1?2α,β)·u(x, t) = f(x, t)，x∈R^d，t∈(0, T)，
u(x,0)=u_0(x)，u_t(x,0)=v_0(x)，lim|x|→∞·u(x, t)=0。
随着分数阶微分在各个领域的应用日益增多，已经相继建立了包含分数阶导数或非局部项的各种偏微分方程(PDE)模型（参见例如[8]、[10]、[20]、[21]、[22]、[23]、[24]）。这些模型的数值解研究进展迅速，但仍有几个突出问题需要进一步关注。已经成功开发了诸如有限差分方法[25]、[26]、有限元方法[27]、谱方法[29]、[30]等技术来处理这些模型。在处理无界域上的PDE时，海mite谱方法特别有利，因为它自然符合无穷远处的衰减边界条件并实现谱精度[31]、[32]、[33]、[34]。得益于这种自然匹配，海mite谱方法已在相关背景下得到广泛应用，特别是在无界均匀介质上的波传播问题中。例如，郭等人[35]构建了一些隐式-显式海mite-Galerkin谱方案，用于模拟具有分数阶拉普拉斯算子的耦合非局部Gordon型系统的动态，并建立了具有稳定性分析的完全离散方案。同时，凌等人[12]、[36]将海mite谱方法应用于扩散粘性波动方程，并为得到的半离散方案提供了误差估计。
尽管模型(4)已经建立多年，但在时间域内直接求解它的方法仍未得到彻底研究。对这一问题的有限数值研究可能源于两个主要困难。第一个困难来自模型中存在的非局部时间算子的数值处理，第二个困难则源于缺乏适当的数学工具或严格分析的理论框架。为了更好地阐明模型的物理意义并扩大其适用性，本文提出了一种在海时间域内直接求解它的海mite谱方法。基于温和分数阶微分，我们将模型(4)重新表述为等效形式(7)。这种重新表述使得可以进行严格的稳定性分析并开发专门的数值方案。在空间上，采用了海mite谱方法。时间步进方法的构建至关重要地依赖于时间项的重新表述：我们首先将方程中的卷积核转换为温和的Caputo分数阶导数的标准形式。然后使用[37]、[38]、[39]中开发的加权移位Grünwald差分算子对当前分数阶导数进行离散化，而二阶时间导数通过中心差分方案近似。最后，提供了方案稳定性和收敛性的严格理论分析。
本文的其余部分结构如下。第2节首先阐述问题的弱形式，并使用海mite谱方法构建空间半离散方案。严格建立了初值问题的弱解的适定性，随后是对半离散方案的详细稳定性分析和误差估计。第3节将讨论扩展到完全离散方案。第4节提供了数值示例，以验证理论结果并展示所提出方案的有效性。最后，第5节总结了本文。

半离散海mite谱形式
我们首先介绍后续分析中使用的函数空间。对于函数u(x)∈L^2(R^d)，其傅里叶变换表示为u^(ξ)。令|ξ|_1、|ξ|_2和|ξ|_∞分别表示ξ的l1、l2和l∞范数。此外，设w(x)是一个权重函数，并定义加权勒贝格空间Lw^2(R^d)={u:R^d→C | ∥u∥_Lw^2<∞}，它包含可测函数，并配备了以下内积和范数（参见[40]）：
(u, v)_{Lw^2} = ∫_R^d u(x)·v(x)·w(x)·dx，∥u∥_Lw^2 = (u, u)_{Lw^2}^1/2。
当w(x) ≡ 1时，我们恢复了完全离散的海mite谱方法

完全离散海mite谱方法
接下来我们讨论完全离散方案。设0=t_0< />< />u^(n+1/2) = 1/2(u^(n+1) + u^n)，?^t_u^n = u^(n+1)?u^n?1/2τ，?^2_t_u^n = 1/τ^2(u^(n+1)?2u^n+u^n?1)。
存在多种离散化(6)中分数阶导数的数值方案。然而，直接近似非局部项2a_0_cdt^(1?α,β)·u(t) + a^2_0_cdt^(1?2α,β)·u(t)是

数值结果
在本节中，通过一个、两个和三个空间维度的综合数值实验证明了所提算法的有效性。下表中的toc表示理论收敛阶数。

示例1
考虑一个涉及单个核函数的一维问题(7)，即
{?^2u(x_1, t) ? ?^1u(x_1, t) + 2a_0_cdt^(1?α,β)·u(x_1, t) = f(x_1, t)，x_1∈r，t∈(0,1)，u(x_1,0)=0，?^u(x_1,0)=0，lim|x_1|→∞·u(x_1, t)=0。
我们取精确解为u(x_1, t)=e^(-βt^4)e^(-x_1^2)，以及相应的

结论
在这项工作中，开发了一种用于模拟微异质多孔介质中波传播的高效海mite谱方法。我们没有直接处理原始模型，而是使用温和分数阶微分得到一个等效模型。所提出的方法制定了半离散和完全离散方案，利用海mite谱方法进行空间离散化，wsg差分公式进行时间分数阶导数的离散化，以及中心差分

未引用的参考文献
[48] credi

作者贡献声明
侯泽鑫：写作——原始草案，可视化，方法论，概念化，软件。
李灿：监督，写作——审阅与编辑，方法论，验证，形式分析，概念化，软件。
郭世民：写作——审阅与编辑，方法论，验证，形式分析，概念化。

利益冲突声明
本手稿的提交不存在利益冲突，所有作者均同意发表。 t)的划分，时间步长为τ=t_n+1?t_n，我们用以下符号表示离散时间级上的函数 u^(n+1 2)=1/2(u^(n+1) + u^n)，?^t_u^n=u^(n+1)?u^n?1/2τ，?^2_t_u^n = 1 τ^2(u^(n+1)?2u^n+u^n?1)。 存在多种离散化(6)中分数阶导数的数值方案。然而，直接近似非局部项2a_0_cdt^(1?α,β)·u(t) + a^2_0_cdt^(1?2α,β)·u(t)是 数值结果 在本节中，通过一个、两个和三个空间维度的综合数值实验证明了所提算法的有效性。下表中的toc表示理论收敛阶数。 示例1 考虑一个涉及单个核函数的一维问题(7)，即 {?^2u(x_1, t) ? ?^1u(x_1, t) + 2a_0_cdt^(1?α,β)·u(x_1, t)=f(x_1, t)，x_1∈r，t∈(0,1)，u(x_1,0)=0，?^u(x_1,0)=0，lim|x_1|→∞·u(x_1, t)=0。 我们取精确解为u(x_1, t)=e^(-βt^4)e^(-x_1^2)，以及相应的 结论 在这项工作中，开发了一种用于模拟微异质多孔介质中波传播的高效海mite谱方法。我们没有直接处理原始模型，而是使用温和分数阶微分得到一个等效模型。所提出的方法制定了半离散和完全离散方案，利用海mite谱方法进行空间离散化，wsg差分公式进行时间分数阶导数的离散化，以及中心差分 未引用的参考文献 [48] credi 作者贡献声明 侯泽鑫：写作——原始草案，可视化，方法论，概念化，软件。 李灿：监督，写作——审阅与编辑，方法论，验证，形式分析，概念化，软件。 郭世民：写作——审阅与编辑，方法论，验证，形式分析，概念化。 利益冲突声明>u^(n+1/2) = 1/2(u^(n+1) + u^n)，?^t_u^n = u^(n+1)?u^n?1/2τ，?^2_t_u^n = 1/τ^2(u^(n+1)?2u^n+u^n?1)。
存在多种离散化(6)中分数阶导数的数值方案。然而，直接近似非局部项2a_0_cdt^(1?α,β)·u(t) + a^2_0_cdt^(1?2α,β)·u(t)是

数值结果
在本节中，通过一个、两个和三个空间维度的综合数值实验证明了所提算法的有效性。下表中的toc表示理论收敛阶数。

示例1
考虑一个涉及单个核函数的一维问题(7)，即
{?^2u(x_1, t) ? ?^1u(x_1, t) + 2a_0_cdt^(1?α,β)·u(x_1, t) = f(x_1, t)，x_1∈r，t∈(0,1)，u(x_1,0)=0，?^u(x_1,0)=0，lim|x_1|→∞·u(x_1, t)=0。
我们取精确解为u(x_1, t)=e^(-βt^4)e^(-x_1^2)，以及相应的

结论
在这项工作中，开发了一种用于模拟微异质多孔介质中波传播的高效海mite谱方法。我们没有直接处理原始模型，而是使用温和分数阶微分得到一个等效模型。所提出的方法制定了半离散和完全离散方案，利用海mite谱方法进行空间离散化，wsg差分公式进行时间分数阶导数的离散化，以及中心差分

未引用的参考文献
[48] credi

作者贡献声明
侯泽鑫：写作——原始草案，可视化，方法论，概念化，软件。
李灿：监督，写作——审阅与编辑，方法论，验证，形式分析，概念化，软件。
郭世民：写作——审阅与编辑，方法论，验证，形式分析，概念化。

利益冲突声明
本手稿的提交不存在利益冲突，所有作者均同意发表。>