基于时间有限元法的矩阵指数鲁棒计算方法

《Mathematics》：A Robust Method for Calculating the Matrix Exponential Based on the Time Finite Element Method

　　

基于常微分方程求解器计算常数矩阵矩阵指数的方案通常被认为不如缩放与平方法（scaling and squaring）等其他策略高效，原因在于前者未能利用待处理矩阵为常数矩阵这一特性。本研究中，研究人员设计了一种基于时间有限元法（Time Finite Element Method, TFEM）的策略，该策略充分利用了上述特性，从而形成一种极为高效的方案——仅需在算法初始阶段一次性求解线性方程组。此外，由于待求逆矩阵即使在被处理矩阵呈现病态条件时仍保持良好条件性（well-conditioned），该策略展现出极高的可靠性，这一特性已通过多个具有挑战性的实例加以验证。本文所述方法可便捷推广至矩阵为时间依赖情形以及存在时间依赖强迫函数的情形。

研究背景与问题阐述

矩阵指数（matrix exponential）在物理系统数学建模中具有核心地位，其计算问题长期困扰数值分析领域。给定常系数线性常微分方程组 ^˙x(t) = Ax(t) 在初始条件 x(0) = x₀ 下的解析解为 e^Atx₀，其中 A 为 n×n 常数矩阵。尽管矩阵指数可通过无穷级数定义直接展开，但截断级数计算往往导致显著误差，数值可靠性极差。Moler 和 van Loan 在其经典综述中系统讨论了缩放与平方法、Chebyshev 有理逼近、常微分方程方法及多项式方法等多种策略，并明确指出：基于常微分方程求解器的方法因未能利用微分方程的线性常系数特性而计算代价高昂——这些求解器反复执行矩阵与向量的乘法运算，因其假设矩阵可能已在各次运算间发生变化。现有指数积分器（exponential integrators）虽有所发展，但其核心仍依赖缩放与平方算法近似矩阵指数。此外，传统方法在矩阵 A 病态或缺乏线性无关特征向量组时频繁失效，凸显该领域对鲁棒算法的迫切需求。

鉴于上述瓶颈，研究人员致力于开发一种既能充分利用常系数特性以提升效率、又具备高度数值稳定性的新型算法，同时要求该算法易于推广至时变矩阵及含强迫函数的非齐次情形。

主要技术方法

本研究的核心技术方法基于时间有限元法（Time Finite Element Method, TFEM），这是一种时域加权残值法，采用分段多项式逼近待求变量，权函数较试探函数低一阶。研究人员分别构建了二次、三次和四次时间有限元格式。关键创新在于：针对常数矩阵 A，通过等时间步长离散，将矩阵求逆运算集中于算法初始阶段一次性完成；后续各时间步仅需执行预计算矩阵与当前时刻指数矩阵的乘法运算。该方法具有能量守恒特性，即使对于保守系统也能在离散层面保持精确守恒，从而赋予算法无条件稳定性。对于二次时间有限元，推导出递推关系 M₂ = [2I ? (Δt/2)A]^?1[2I + (Δt/2)A]，其中 I 为单位矩阵；三次和四次格式分别对应 M₃ 和 M₄ 的类似表达式。算法实现中，仅需在初始时刻求解含 n 个右端项的线性系统以获得 M，随后通过递推 E_j+1 = MEj 逐步推进，直至最终时刻 T。该方法的显著优势在于：即便 A 病态或接近奇异，由于单位矩阵 I 的存在，M 相关矩阵始终保持良好条件性；且算法不依赖 A 的任何特殊结构（如对角化、Jordan 标准型或特征向量完备性），适用于任意矩阵。

研究结果介绍

截断级数失效验证：针对文献中截断无穷级数导致误差的经典反例，采用 Δt = 1/16 的四次时间有限元计算，结果与文献值及 Mathematica 参考解在 8 位小数上完全一致，误差范数随步长减小呈现系统性收敛。

与解析解对比：对具有已知解析解的矩阵，采用 Δt = 1/16 计算，所提算法结果与解析解几乎完全吻合，仅部分元素在第 14 至 15 位小数出现偏差。

近亏损矩阵处理：对于特征向量方法易失效的近亏损（nearly defective）矩阵，取 Δt = 0.1 时误差范数达 O(10^?10) 量级，展现出对较大步长的高精度容忍度。类似地，针对文献报道的另一近亏损实例同样获得正确解。

Jordan 块矩阵求解：针对 Jordan 块形式的特殊矩阵（其矩阵指数具有已知解析表达式），取 Δt = 0.1 时误差范数约为 O(10^?12），近乎完美复现理论解。

病态矩阵挑战：针对高度病态矩阵（旧版 MATLAB expm() 曾给出错误结果的实例），取 Δt = 0.1 计算所得结果与当前版本 MATLAB 一致，误差范数随步长减小而改善。另一样例中矩阵元素变化幅度极大，所提算法结果与 MATLAB expm() 高度吻合，而基于特征值/特征向量的 expmdemo3 函数因病态问题给出警告及错误结果。

反对称矩阵正交性保持：针对已知指数矩阵应为正交矩阵的反对称矩阵实例，观察到需极小子时间步长以捕捉解的高频振荡特性。此例中尺度与平方法（expm()）在精度和计算时间上优于所提算法，这归因于时间步进方法的误差累积效应。

大尺度奇异矩阵计算：针对 100×100 全 1 矩阵（秩 1，仅一非零特征值 100）这一大规模奇异实例，Fortran 实现的缩放与平方法耗时约 0.03 秒得对角元约 2.6949×10⁴³，而所提算法取 Δt = 0.001 和 Δt = 0.0005 分别耗时 0.70 秒和 1.10 秒，获得精确得多的结果 2.7183×10⁴³。若需获取 1000 个时间步的完整动力学解，所提算法的逐步推进特性反而具有整体效率优势。

三角矩阵函数计算：以实例 2 的矩阵为基础，分别计算矩阵余弦和矩阵正弦，取 Δt = 1/16 时与解析解几乎完全匹配，验证了方法对一般矩阵函数的适用性。

讨论与结论总结

本研究开发了一种鲁棒的矩阵指数计算策略，其核心在于通过时间有限元递推关系实现逐步计算。该算法无需用户调整任何内部参数或外部因子，具有显著的实现简洁性。通过与现有策略的性能对比，该方法在一系列具有病理特征（pathological）的测试实例上展现出卓越的数值鲁棒性。

相较于传统技术的核心优势体现在两个维度：其一，算法能够计算矩阵 A 在完整时间区间 [0, T] 上的全部时间演化过程，而非仅获取单一时刻的指数矩阵；其二，与传统方法形成鲜明对照，所提方法可便捷推广至 A 为时间函数且控制方程含强迫函数的非齐次情形（即 ^˙x(t) = A(t)x(t) + f(t)），这一扩展能力在现有主流方法（如缩放与平方法）中难以实现。论文发表于《Mathematics》期刊。