一种基于Crank-Nicolson ADI Galerkin方法的数值模拟技术，用于求解粘弹性板中出现的非线性Volterra积分微分模型

《Applied Mathematics and Computation》：A Crank-Nicolson ADI Galerkin approach for numerical simulation of nonlinear Volterra integrodifferential models arising in viscoelastic plates

　　胡法江|埃姆兰·托希迪|达旭
湖南工业大学理学院，中国湖南省株洲市412007

摘要
本文提出了一种高效且强大的方法，用于处理非线性Volterra积分微分方程，这些方程在研究由各向同性粘弹性材料构成的板状结构时尤为常见。我们重点关注那些没有显式解析表达式的卷积核，而是通过它们的拉普拉斯变换来描述。所提出的方案结合了Crank-Nicolson方法和二阶卷积积分技术以实现时间上的局部化，并在空间变量上采用有限元Galerkin逼近。为了降低多维问题中的计算成本，引入了ADI（交替方向隐式）方法。通过离散能量分析证明了该方法的收敛性和稳定性，在适当的假设下，误差估计为O(τ^2 + hr + 1)，其中h和τ分别表示空间和时间网格的大小。数值实验结果验证了该方法的鲁棒性，并确认了预测的收敛速率。

引言
受粘弹性各向同性板状结构线性理论框架应用的启发[1], [2]，我们研究了以下有趣的数学模型（该模型包含关于时间变量的偏导数和积分项）：
?u(t,χ) ? Δu(t,χ) ? ∫_0^t K(t?s) Δu(s,χ) ds = g(u(t,χ))，
(t,χ) ∈ (0,T] × Ω，u(0,χ) = u_0(χ)，x ∈ Ω，u(t,χ) = 0，(t,χ) ∈ [0,T] × ?Ω。
这里χ = (χ1,χ2)，?tmu := ?mu/?tm（m=1,2,3），Ω = I×I（I=(0,L)），初始数据u_0和非线性激励函数g(u)属于希尔伯特空间H。函数K是从R^+=[0,∞)到R的映射，并且根据[3]假设为局部绝对连续的。对于任何函数Θ，其拉普拉斯变换表示为L{Θ}或等价地Θ^：Θ^(z) = ∫_0^∞ e^(-zt)Θ(t) dt，Θ^(iτ_0) = lim_z→iτ_0, |?(z)|>0 Θ^(z)，τ_0∈R, |?(z)|>0，当极限存在时定义（?表示复数的实部）。从与之相关的拉普拉斯变换K^的角度来看，可以进一步阐述关于K的假设。具体来说，我们的分析基于以下假设[3], [4]：
K^(η) = f(η)a^(η)，|?(η)|>0，
其中f(η) = F((a^(η))^(-1)b^(η))，F(ξ) = m + ξ^q + pξ^q，p>0，m>0，
并且b(t)以及a(t)是凸函数、非负的、连续的、递减的，并且在[0,∞)上不是常数；
此外还有dd′(b(t)) + a′(t))|τ=0 > ?∞的假设。同时假设：?db(τ)/dt和?da(τ)/dt在R^+上具有凸性，非线性激励函数g也应该具有连续的Lipschitz性质。换句话说：∥g(t,x, w_1) ? g(t,x, w_2)∥ ≤ L∥w_1 ? w_2∥，?L>0，w_1,w_2∈L^2(Ω)。

K(t)在(1)中的假设是由于该方程广泛适用于具有粘弹性特性的物质中的热传输等现象[1], [2], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14]而提出的。这些应用特别包括粘弹性各向同性杆的线性理论框架，这对应于特定情况g(u(t),t) = g(t)[15], [16], [17], [18], [19], [20], [21]。此外，方程的核函数是由压缩和剪切的应力松弛模量（a^(η), b^(η)）以复杂的方式组合而成的[3]，其中f(η)表征了压缩和剪切响应之间的频率依赖耦合比。我们在(84)和(87)中给出了它们的具体表达式，从而得到了核K(t)的相应拉普拉斯变换。通过将其与弹性对应物的“对应原理”相结合，可以得到控制线性粘弹性动态问题的方程[15], [20], [22], [23]。我们首先对方程(1)进行形式上的微分（Δu(x,s) = Δ2u(x,s)），然后使用(3)中的参数应用傅里叶变换，得到粘弹性杆振动的转换方程（例如参见[15]，第109-112页）。类似地，与(2)和(3)相关的杆中波的公式也进行了研究，其中参数q=1, p=2, m=0。对于纵波，我们使用Δu(x, s)；而对于弯曲波，则需要Δ2u(x, s)。在所有讨论的情况下，都需要应用自伴类型的边界条件。此外，模量b(t)和a(t)通常表示为衰减指数项的非负线性组合[3]。这种表示方法，结合假设(2)–(6)提供的K函数选择，允许分析方程(1)中解的规则性及其数值近似的误差估计。

现有研究强调了使用卷积函数K的拉普拉斯变换K^(z)的有效性。Xu[4]通过在拉普拉斯域而不是直接使用核函数本身，证明了线性问题(1)的数值近似的均匀?1-收敛性。随后，Qiu等人[24]通过K^(z)分析了K及其导数的性质，显示了它们在数值分析中的关键作用，尽管他们的结果仅限于半离散设置。另一方面，交替方向隐式（ADI）方法由于其解耦多维算子和降低计算复杂性的效率而在高维问题中证明是不可或缺的[25], [26], [27], [28], [29], [30], [31]。受到这些发展的启发，我们旨在弥合半离散理论与实际多维计算之间的差距。具体来说，我们设计了一种完整的离散化Crank-Nicolson (CN) ADI Galerkin方法，将基于拉普拉斯变换的核分析与ADI分割的效率结合起来，并通过能量技术证明了其稳定性和收敛性。

我们的讨论按照以下大纲展开：第2节介绍了构建时间和空间变量离散化的几个前提和基本概念。第3节发展了用于求解所考虑模型的CN ADI Galerkin数值框架。第4节提供了所提方案稳定性的分析和误差估计。最后，第5节提供了验证我们理论结果的测试问题。

几个前提
本节提供了实现问题(1)的空间和时间离散化的几个前提和基本概念。首先，我们将时间变量[0, T]区间离散化为网格点0=t_0 < t_1 < ? < t_N = T，其中tm=mτ，0≤m≤N。假设tm?1/2=tm+tm?1/2，δt∈(0,1)，V_m=V_{m?1}?V_{m?1}τ，1≤m≤N。此外，可以构造L^2(Ω)中的内积及相关范数：
(V,W)=∫_Ω V(x_1,x_2)W(x_1,x_2) dx_1 dx_2，∥V∥=(V,V)。
假设Hm(Ω)是Sobolev空间。

Crank-Nicolson ADI Galerkin方法的表述
为了解决(1)，本文介绍了完整的离散化Crank-Nicolson ADI Galerkin公式。

离散化方法的分析
在本文的这一部分，我们研究了Crank-Nicolson ADI Galerkin方法的稳定性并推导了误差估计。为此，我们首先将方案表示为以下等价形式：
(δt/1)U_1,V_h) + (?U_1?V_h) + ω_0(τ)(?U_1?V_h) = 2ω_0(τ) ? ω_1(τ) + (4/?U_0?V_h) + 1/2(g_1(U_0) + g_0(U_0), V_h)，
(δt/1)U_n,V_h) + (?U_n?1/2,?V_h) + ∑_(p=1)^n ω_n^(-p)(?U_p?1/2,?V_h) = ω_n^(-1)(τ) ? ω_n(τ) + (4/?U_n^(-1)(U_n?1),V_h) + 1/2(g_n(2U_n?1?U_n?2),V_h)，n=2,3,…,N。
然后，定义V_{n}^′=e^(-γtn)V_n，1≤γ<∞，ω_{n}^*(τ)=ω_n(τ)e^(-γt)。

针对测试问题的建议方法实施
所有数值模拟都是使用由方程(25)和(26)定义的CN ADI Galerkin方法进行的，以近似问题(1)的解。代码在MATLAB R2021b中实现，并在配备Intel Core i5-12500H处理器（3.10 GHz）和16.0 GB RAM的Windows 10（64位）台式计算机上执行。所提出的测试问题用于验证方案的效率并评估其数值精度。我们假设I=[0,1]且M_x1=M_x2=M。ADI方法也...

结论
我们通过三个部分总结我们的工作：
(i) 主要结果：在这项工作中，我们提出并分析了一种高阶方案，用于处理由隐式定义核的粘弹性模型引起的Volterra积分微分方程。该方法结合了Crank–Nicolson方案、二阶卷积积分和有限元Galerkin逼近，以及提高多维效率的ADI策略。离散能量分析证明了稳定性以及最优误差估计O(τ^2 + hr + 1)。

利益冲突声明
作者在此声明不存在利益冲突。