基于滤波级数的分数阶数值近似方法及其在计算力学中的应用

《Communications in Nonlinear Science and Numerical Simulation》：Fractional-Order Numerical Approximation Method Based on Filter Series and Its Applications in Computational Mechanics

　　

摘要

分数阶微积分作为整数阶导数的扩展，具有整数阶导数所不具备的时间维度非局部性和记忆特性。它被应用于描述生物组织、矿物和复合材料的粘弹性力学响应，以及控制系统的设计。卷积方法通常用于分数阶微积分的数值实现，但通过该方法实现的固体本构模型存在记忆开销随模拟时间增加的问题。分数阶微积分的另一实现路径是在拉普拉斯域中采用有理分式逼近，由于其对控制系统设计的指导意义而广受欢迎。本文基于有理分式逼近提出了分数阶微积分的时域数值格式。与卷积方法相比，该方法在相同数值精度下显著减少了时间和记忆开销。本文进一步讨论了生物材料中分数阶本构模型的起源，它源于材料的无限阶自相似结构。与几何分形不同，这种自相似特征体现在粘弹性响应算子中，被称为物理分形。利用有理逼近方法，本文在时域中数值实现了具有分数阶特征的致密骨粘弹性本构模型，为进一步的有限元模拟奠定了基础。

引言

作为整数阶微积分理论的自然延伸，分数阶微积分因其对非局部依赖性和历史记忆效应的精确表征能力，已成为描述复杂力学系统（如粘弹性材料和生物组织）行为的核心理工具。它在粘弹性材料力学行为分析、生物神经网络动力学建模和结构动力学响应预测等交叉学科方向展现了独特优势[1]–[10]。然而，分数算子在时间维度的非局部响应特性给数值实现带来了困难，分数算子的高效数值逼近技术已成为该领域从理论建模到工程应用的核心瓶颈[11]。围绕这一关键问题，研究者们提出了各种逼近方法[12]–[16]。

以Grünwald-Letnikov卷积为代表，该方法从分数阶微积分的积分形式出发，将其实现为数值卷积，广泛用于分数阶微积分的数值实现。Grünwald-Letnikov卷积格式是一种常见的离散格式：

$$f(nh) = \sum_{k=0}^{n} w_k f(nh - kh); \quad w_0 = 1, \quad w_k = \left(1 - \frac{1+\alpha}{k}\right) w_{k-1}, \quad \forall k = 1, 2, \cdots$$

然而，此类方法通常具有线性记忆复杂度，当应用于大规模数值计算时，记忆开销随模拟时间增加[17]。

另一种方法是在频域中对分数算子的核函数采用有理分式逼近，由于其对控制系统设计和分数阶电路（AFCs）模拟的指导意义而广受欢迎[18]–[23]。早期逼近方法大多依赖于$p^\alpha$函数的某些级数展开。例如，Padé逼近方法通过构造有理函数实现非线性逼近。尽管计算过程简洁，但精度高度依赖于频段范围。一旦试图拓宽有效频段，由于零极点分布未优化，其误差（尤其是相位误差）将急剧（指数级）增加，难以满足宽频域力学模拟的需求[24]；连分数展开法将分数算子转化为递归的分数结构以实现逼近阶数的灵活调节，但截断误差与应用频段宽度显著正相关。尽管高阶截断可在一定程度上拓宽有效频段，但会导致系统结构复杂度呈几何级数增长，显著增加了数值模拟和工程实现的难度[25]。Carlson[26]提出了一种迭代方法来逼近分数算子。该方法获得的有理分式在中低频段具有良好的逼近效果和结构稳定性，但其极点分布中部密集而两端稀疏。高频段和极低频段的稀疏极点导致逼近精度迅速衰减，难以满足宽频域电路模拟的需求。Oustaloup提出了一种可人工设计工作频段的有理分式逼近。通过在预设目标频段内对零极点进行对数均匀排列，该方法实现了分数算子的高精度逼近。其核心优势在于强频段可调性，可根据实际应用场景需求灵活设置截止频率范围，并在宽频域内保持稳定的幅相特性，为宽频域分数阶电路系统模拟提供了可能[27]；Matsuda逼近法采用折线拟合策略来逼近理想分数算子的Bode图特性，通过对幅相特性曲线进行分段线性化简化计算。该方法计算复杂度和工程实现成本低，但逼近精度受限于分段数。高频段容易发生相位失真，影响分数阶系统的动态响应模拟精度[28]。系统性性能对比分析表明，Oustaloup逼近法在频段可调性和宽频域逼近精度方面表现最佳，更适合电路设计和计算力学等领域宽频域高精度分数阶系统的建模需求[11], [29]–[31]。其独特的拓扑结构和极点配置原则为本文提出的滤波算子化思想奠定了重要基础。

生物材料力学旨在连接力学性能、生物功能、几何结构和组成成分，揭示生物材料卓越性能背后的力学和物理机制[32]。生物组织和结构从软质的皮肤、肌肉、毛发到坚硬的骨骼、鳞片、贝壳和外骨骼系统，通常具有适应其生理功能的分层复合架构[33]–[35]。典型的天然复合材料，如珍珠层、肌腱和蜘蛛丝，通过软硬自组装一体化策略实现有利的力学性能和生物功能[36]–[38]。这些系统的共同特征是高程度自相似性。早期研究通常使用欧几里得空间中的经典几何分形来描述这一特征[39]–[41]。然而，对生物粘弹性的纯几何分形描述仍然有限，因为它仍然依赖于主要为均质无机材料开发的连续介质假设，倾向于将问题简化为实验拟合和现象学描述。因此，它可能忽略了多尺度结构与遗传响应之间的内在联系。尽管准线性粘弹性和多元素Kelvin-Voigt模型或其非线性扩展可以拟合特定的实验曲线[42], [43]，但这些模型本身并不能解释材料自相似性如何产生分布式松弛。同样，仅将分数阶导数作为拟合工具使用或不与材料层级建立物理联系的拓扑模型构建[44]，也无法完全克服结构与功能之间的分离。

近年来，为了解决上述问题，研究者们提出了"物理分形"理论。以真实生物体的自相似结构为中心，该理论为深入理解生物材料等复杂力学系统的粘弹性行为提供了新的视角[45]–[50]。人体许多组织和结构具有天然分形特征（如致密骨纤维[45]–[49]、肌肉/韧带纤维、神经纤维[50]以及动脉血流等）。这些天然分形结构导致其粘弹性响应表现出显著的分数阶特征，难以通过传统整数阶模型精确表征，而分数阶力学系统能够完美匹配这一特征。基于物理分形的思想，我们提出通过连分数展开法将分数阶算子转化为弹簧-阻尼器网络的递归拓扑结构[45]–[47]，实现生物力学行为的数值模拟和物理逼近，为生物医学仿生等计算力学应用领域提供支持。然而，连分数截断方法[48]在计算力学的实际应用中存在固有局限性：首先，截断连分数算子的极点分布由数学展开形式固有决定，无法根据目标力学信号的频段特性灵活调节。当输入信号处于极高或极低频段时，稀疏的极点导致逼近效果不佳，无法精确匹配生物材料等复杂力学系统的实际响应特性；其次，连分数截断级数与有效工作频段之间的相关性缺乏明确的定量规律。为覆盖所需工作频段，需要盲目增加截断级数，导致系统部件数量和计算功耗呈指数增长，显著增加了力学系统数值模拟的成本；第三，连分数展开得到的算子参数为固定值，缺乏动态调节能力，难以适应生物材料在不同生理状态下差异化的频段响应需求。这些局限性严重制约了分数阶力学系统在生物医学仿生和宽频域结构模拟等高端计算力学领域的应用拓展。在此背景下，基于算子化思想的滤波算子[49]为解决这一计算力学核心瓶颈问题提供了有前景的技术方案。

首先，本文基于有理分式逼近提出了分数阶算子的时域数值实现方法。与卷积方法相比，该方法显著节省了计算时间和记忆开销。对于一类具有半阶分数阶特征的算子，本文提出了一种基于将此类算子表示为奇异积分的新有理逼近方法，优于现有的Oustaloup有理逼近。在本文后半部分，我们讨论了分数阶本构模型的物理起源，它源于材料的无限阶自相似结构。与空间分形不同，这种自相似特征被抽象为结构的粘弹性响应算子，形成算子域中的分形，因此被称为分形算子。本文进一步讨论了使用有理逼近方法实现该粘弹性本构模型的可能性，并给出了致密骨的数值实现。需要注意的是，本文采用的皮质骨模型作为概念性和示范性示例。其目的是说明基于滤波序列的有理逼近方法构建分数阶粘弹性本构模型的可行性，而不是提供皮质骨实际力学行为的全面定量表征。本工作的核心贡献在于分数阶算子的通用数值逼近框架，而骨模型仅作为该框架的典型应用案例呈现。

章节片段

分数阶算子有理分式逼近

从Liouville形式的分数阶算子定义出发，连续卷积可在时间维度上离散化以实现分数阶算子的数值实现。Grünwald-Letnikov卷积格式(1)通常作为分数阶本构模型中分数阶算子数值实现的基础。该方法需要存储整个激发历史，因此空间开销为O(L)（其中L为时间步数）。

算子化表示

基于对分数阶算子的近期研究，我们提出了一种基于算子化思想的分数阶系统设计框架[49]，其中包括算子域的重定义。针对Mikusinski单一算子域概念的模糊性，我们重建了运算微积分的代数基础，建立了以微分算子p为核心的算子域和以卷积为核心的卷积域平行对偶的双域结构。

讨论

所提出的滤波序列实现方法应被视为分数阶或分布式松弛核的计算表示。其主要优势在于记忆项通过局部内部变量进行更新，而非存储整个应变历史。这使得该方法对于有限元模拟具有吸引力，因为在有限元模拟中，分数阶本构更新必须在每个积分点进行。

每个积分点的内存开销为 O(N)，每个时间步的计算量为 O(N)。

**局限性**与**未来工作**

首先，当前的数值演示侧重于算子逼近和计算效率；它们尚未构成针对特定骨皮质骨实验的完整验证。其次，Oustaloup 及相关有理逼近方法需要预选频率范围；在该范围之外响应可能恶化，因此输入频谱估计和收敛性检查至关重要。第三，公式 (8) 中使用的线性插值需要足够小的时间步长，以便……

**结论**

本文提出了一种基于算子域中有理分式逼近方法的高效分数阶微积分数值实现方案，在时间成本和内存开销方面显著优于卷积方法。

有理逼近方法在频域中的逼近效果特征为：在有效频带内性能优异，而在频带之外则失效。本文认为，该有效频带……

**声明**

- **资助来源**：本工作得到了国家自然科学基金（项目编号：12050001）和成都理工大学人才项目（项目编号：2025RC033）的资助。
- **利益冲突**：不适用。
- **伦理审查及参与同意**：不适用。
- **出版同意**：所有作者均已阅读并批准最终手稿。
- **数据可用性**：可向通讯作者索取。
- **材料可用性**：可向……索取。

**CRediT 作者贡献声明**

**笪智默**：写作 – 审阅与编辑，写作 – 原创草稿，验证，方法论，调查研究，形式分析，概念化。
**李浩**：调查研究，形式分析，概念化。
**蒋智真**：软件，方法论，形式分析。
**彭刚**：写作 – 审阅与编辑，调查研究，形式分析，概念化。
**蒋瑞恒**：调查研究。
**周天一**：调查研究。
**尹亚军**：写作 – 审阅与编辑，写作 – 原创草稿，方法论，概念化。

**利益冲突声明**

作者声明，他们不存在任何已知的竞争性财务利益或个人关系，这些利益或关系可能会影响本文所报告的工作。

笪智默 | 李浩 | 蒋智真 | 彭刚 | 蒋瑞恒 | 周天一 | 尹亚军

清华大学工程力学系，北京，100084，中国