全局时间分数阶扩散的快速正弦变换预条件方法

《Fractal and Fractional》：Fast Sine-Transform Preconditioning for Global-in-Time Fractional Diffusion

　　

时间分数阶扩散方程描述了异构介质中的亚扩散输运过程，但由于非局部Caputo导数以及解在初始时刻产生的弱奇异性，其数值处理变得复杂。研究人员研究了一种全局时间离散化方法，该方法将时间域的谱配置与空间域的二阶守恒有限差分模板相结合。其中，时间域采用在Chebyshev–Gauss–Lobatto节点上求值的分数次幂基{t_l^α}_l=0^N，该基能够再现解的奇异展开的主导项；空间域则在单元面上使用扩散率的调和平均，从而在间断介质中保持精度。由此产生的全离散问题是一个大型、非对称且在时间上稠密的线性系统，其二范数条件数随空间步长的平方倒数增长，因此无预条件的Krylov子空间迭代在网格细化时会发生停滞。利用离散算子的Kronecker和结构，研究人员通过离散正弦变换 (DST) 对空间因子进行快速对角化来构建预条件子。对于常数扩散率，该预条件子能精确再现算子并产生一个直接求解器；对于变扩散率，它在谱上等价于算子，研究人员证明了预条件系统的特征值聚集在一个以1为中心的圆盘内，其半径仅取决于系数对比度，而与网格、时间自由度数或分数阶无关。一维和二维空间数值实验证实了二阶空间精度，且预条件迭代次数保持平坦——一维中从M=32到M=8192保持12次迭代，二维中每个方向最高达到M=512时保持11次——而无预条件的迭代次数则增长超过两个数量级。在时间域，精度呈谱精度，直到幂基病态范德蒙德 矩阵的舍入误差占据主导：该限制在特定自由度N和分数阶α下达到，此时可达到的误差约为10^-6。通过与均匀网格和分级网格上的L1格式、Alikhanov L2-格式以及Grünwald–Letnikov卷积求积进行基准测试，量化了全局方法的适用时机：在强制问题和具有适中刚度的模式上，该方法达到规定精度的速度快一到两个数量级，且内存消耗减少数倍；而对于强衰减模式，分数次幂基仅呈代数收敛，在相对误差低于某特定值时，分级时间步进更为可取。

在物理与生物系统的诸多领域，如多孔或裂隙岩石中的输运、非晶半导体中的电荷传导以及拥挤细胞质中的蛋白质运动等，常观察到反常扩散现象，其均方位移随时间呈次线性增长。通过使用Caputo分数阶导数替换经典热传导方程中的一阶时间导数，可以建立时间分数阶扩散方程模型。该方程在时间上具有非局部性，状态通过弱奇异记忆核依赖于过去的全部历史，这种记忆特性使得模型忠实于亚扩散动力学，但也导致求解过程计算代价高昂。

目前数值处理面临两大难题：首先，即使输入数据光滑，解在初始时刻t=0处也不光滑，存在弱奇异性。这导致基于多项式插值的传统时间步进格式（如L1格式）在原点附近精度下降，除非采用非均匀的分级网格或进行奇异项修正；其次，非局部性使得每一时间步均与之前所有历史步耦合，传统步进法的计算成本随步数呈二次方增长。若采用全局时间离散化，虽能提供高精度，但会产生大型、非对称且在时间上稠密的线性系统，其条件数随空间网格的细化而急剧恶化，导致无预条件的Krylov子空间方法在网格细化时收敛停滞。针对这些问题，研究人员开展了一项研究，旨在通过挖掘空间算子的可分离结构来构建高效预条件子，从而解决全局配置法的计算瓶颈。该研究得出结论，利用离散正弦变换构建的预条件子能够使Krylov迭代次数完全独立于空间网格和时间自由度，这项研究发表在《Fractal and Fractional》上，为求解分数阶扩散方程提供了一种精度高、成本低且结构易于并行化的实用方法，具有重要的学术与工程意义。

研究人员为开展此项研究采用的关键技术方法包括：在时间域采用基于Chebyshev–Gauss–Lobatto节点的分数次幂基{t_l^α}_l=0^N进行谱配置，以精确捕捉初始弱奇点；在空间域采用基于调和平均扩散率的二阶守恒有限差分模板；将全离散问题组装成具有Kronecker和结构的线性系统。在求解器方面，采用广义最小残差法 (GMRES) 结合基于离散正弦变换 (DST-I) 的快速对角化预条件子进行迭代求解。研究基于纯数值模型构造算例，未涉及具体生物样本队列。

精度与标准时间步进格式的比较
研究人员通过构造半离散制造解来研究精度。时间域上，由于分数次幂基能够吸收初始奇点，方法展现出谱精度，直到病态范德蒙德 矩阵的舍入误差占据主导。研究确定了该精度壁垒，指出对于特定分数阶，可用时间自由度存在上限，超出此界限增加自由度反而恶化结果，此时可达到的误差约为10^-6。空间域上，通过选择时间精确的制造解，证实了方法在所有测试分数阶下均保持二阶空间精度。在与均匀网格及分级网格上的L1格式、Alikhanov L2-格式以及Grünwald–Letnikov卷积求积等标准时间步进格式的基准测试中，研究量化了全局方法的适用范围：对于强制问题和中等刚度模式，全局方法在达到规定精度时速度快一到两个数量级，且内存减少数倍；但对于强衰减模式，分数次幂基仅呈代数收敛，此时分级时间步进在相对误差低于特定值时更具优势。

求解器性能与鲁棒性
在空间网格细化实验中，未预条件的GMRES迭代次数随条件数恶化而急剧增长，而采用快速对角化预条件子的GMRES迭代次数保持平坦（一维12次，二维11次），证实了网格独立性。研究还证实了重启参数对预条件迭代无关。鲁棒性实验表明，预条件迭代次数与分数阶无关，且随系数对比度的增加而缓慢增长。研究人员直接验证了特征值聚类定理，预条件算子的特征值确实聚集在以1为中心的圆盘内，且观测到的最大偏差与预测半径高度吻合。在不连续系数实验中，即使跳跃比高达100，调和平均与算术平均的迭代次数也表现一致且收敛迅速，表明该预条件子在极端对比度下仍能优雅退化，这归因于特征值在圆盘内高度聚集而非均匀填充。

在讨论部分，研究人员总结了研究的局限性并指出了未来方向。研究结论部分翻译如下：研究人员研究了一种时间分数阶扩散方程的全局时间离散化方法，该方法将时间域上适用于弱初始奇点分数次幂基的谱配置，与空间域上的二阶守恒有限差分模板相结合。该离散化产生了一个非对称且病态的Kronecker和线性系统，其二范数条件数随空间步长平方的倒数增长。核心贡献是通过离散正弦变换对空间因子进行快速对角化而获得的预条件子：对于常数扩散率，它是精确的直接求解器；对于变扩散率，它是谱等价的，因而是与网格无关的预条件子。研究人员证明了预条件谱聚集在以1为中心的圆盘内，其半径等于系数对比度，与空间网格、时间自由度数和分数阶无关。一维和二维数值实验与此预测相符：一维在空间细化下保持12次平坦迭代，二维未知量最高达约262,000时保持11次，而无预条件次数增长超过两个数量级。与L1、分级L1、L2-和卷积求积的基准测试确定了该离散化的定位：在强制问题和中等刚度模式下，它达到规定精度的速度快一到两个数量级，内存消耗少数倍；对于强衰减模式，分数次幂基仅代数收敛，分级步进更可取。
此外，研究人员指出存在三大限制：一是单项式基的范德蒙德矩阵条件数随规模几何增长，存在硬计算壁垒；二是需引入正交基（如广义Jacobi函数）以根本消除该壁垒；三是分数次幂基无法表示强衰减模式的代数尾部，需采用多区域配置。未来工作中，GPU加速和针对非线性方程的Newton-Krylov循环扩展是极具潜力的方向。