针对高振荡常微分方程的一类新型指数积分器的误差分析——基于局部线性延拓技术

《Mathematics and Computers in Simulation》：Error analysis on a novel class of exponential integrators with local linear extension techniques for highly oscillatory ODEs

　　

# 摘要

本文研究了一类具有有界振荡能量的非自治常微分方程，其线性部分的特征值为纯虚数且大小与小参数?成反比。当0 0$为不依赖于$\varepsilon$的终止时间。常矩阵$\mathbf{A} = (A_{ij})_{d \times d}$假设为可对角化的且特征值为纯虚数，但不一定是对角阵。依赖于时间的非线性强迫项$\mathbf{F}$独立于$\varepsilon$。当$0 < \varepsilon \ll 1$时，系统(1.1)描述了由线性振荡动力学主导的高度振荡问题。一个典型示例是二阶振荡方程[5]–[7]：

$$\ddot{\mathbf{y}} + \frac{1}{\varepsilon^2}\mathbf{M}\mathbf{y} = \mathbf{g}(\mathbf{y}, t), \quad t \in [0, T], \quad \mathbf{y}(0) = \mathbf{y}_{\text{in}}, \quad \dot{\mathbf{y}}(0) = \dot{\mathbf{y}}_{\text{in}},$$

其中$\mathbf{M}$为对称正定矩阵，$\mathbf{g}$为非线性函数。引入缩放动量变量$\mathbf{p} = \varepsilon\dot{\mathbf{y}}$，令$\mathbf{u} = [\mathbf{y}^\top, \mathbf{p}^\top]^\top$，从而将系统(1.2)改写为(1.1)的形式。此类二阶高度振荡系统在科学计算中广泛存在，出现在众多学科中，包括经典力学、分子动力学、量子物理、天文学、电路仿真和等离子体物理[8]–[10]。本文主要研究应用于(1.1)的数值方法，其背景为(1.2)。将相关讨论推广到(1.1)的更一般情形是直接的。

若$\mathbf{g}$不依赖于时间，且存在势函数$U(\mathbf{y})$满足$\mathbf{g}(\mathbf{y}) = -\nabla U(\mathbf{y})$，则(1.2)表示一个Hamilton系统[11]。Hamilton函数取如下形式：

$$H(\mathbf{y}, \dot{\mathbf{y}}) = H_1(\mathbf{y}, \dot{\mathbf{y}}) + U(\mathbf{y})$$

其中

$$H_1(\mathbf{y}, \dot{\mathbf{y}}) = \frac{1}{2}\|\dot{\mathbf{y}}\|^2 + \frac{1}{2\varepsilon^2}\mathbf{y}^\top\mathbf{M}\mathbf{y},$$

$\|\cdot\|$表示标准Euclidean范数。由(1.4)定义的$H_1$代表(1.2)中线性振荡器分量所对应的总动能与势能。

具有有界能量的系统无论在理论研究还是应用研究中都受到广泛关注，特别是在经典力学中[8]、[12]。更一般地，(1.2)的右端项可能同时建模内部力和外部时变相互作用，将(1.2)的适用范围扩展到保守系统之外。尽管如此，振荡能量(1.4)的定义仍然适用。由于在高度振荡问题中$H_1$通常主导总能量$H$，因此我们引入有界振荡能量条件：

$$H_1(\mathbf{u}) \approx \mathcal{O}(1), \quad t \in [0, T],$$

作为更一般的能量约束。该条件表明系统中的振荡幅度较小。

从计算角度来看，系统(1.1)中占主导地位的线性项对常规数值方法（包括标准Runge–Kutta和Runge–Kutta–Nystr?m格式）施加了严苛的稳定性约束[1]、[7]。为保证数值可靠性，时间步长$h$必须分辨最快的振荡，从而需要满足限制条件$h < \mathcal{O}(\varepsilon)$[12]、[13]，这限制了实际计算效率。此外，大多数常规数值方法的误差随频率增大而增大。其构造通常基于精确解的Taylor展开，其余项通常依赖于$\varepsilon$的负幂次[14]、[15]，导致关于$\varepsilon$不具有均匀精度。

对(1.2)的数值方法研究可追溯到20世纪60年代。Gautschi首次引入三角多项式作为滤波器来处理高频振荡[16]。此后，各种Gautschi型三角积分器被相继发展[17]–[19]。后续研究考察了长时动力学意义上的数值解[20]–[22]，以及Hamilton系统背景下的数值格式[13]、[17]、[19]、[20]、[23]、[24]。若干类方法通过二阶方程的常变分公式构造[23]–[25]。这些处理在适当的步长条件下显著改善了精度和稳定性。

近几十年来，指数积分器（EIs）已成为求解一阶刚性微分方程的强大工具[26]–[31]。EIs的基本构造利用了常变分公式。对于刚性主要集中在线性部分的系统，由于其优异的稳定性特性，该方法尤为 advantageous[26]。后续发展还将EIs与线性化技术相结合[29]、[32]–[34]，例如著名的指数Rosenbrock型方法[27]、[35]、[36]。这些方法在每个时间步长上将非线性项的Jacobian矩阵并入线性部分后再进行指数积分。然而，这些传统指数型方法要求刚性系统的解具有时间光滑性[31]，在$\varepsilon \ll 1$的参数区域下该条件通常被违反，导致均匀精度的丧失。

本工作聚焦于提升应用于具有有界振荡能量的非自治系统(1.1)的EIs的线性化技术。所提出的方法解决了传统EIs应用于(1.1)时的关键局限，包括收敛所需的步长限制和缺乏均匀精度。因此，所提格式即使对于大于$\varepsilon$量级的步长仍能保持其收敛阶。此外，新的EIs不需要阶数条件，使我们能够建立构造任意高阶格式的统一框架。这些特性使得该方法对于经典力学中具有有界能量的问题尤为有效。

我们的数值格式通过引入辅助变量来增维原始系统。在过去十年中，辅助变量方法得到了快速发展，例如广为人知的标量辅助变量（SAV）[37]–[39]方法，它以较低计算代价促进了保持能量耗散或守恒的格式构造。相比之下，我们的增维方法动机源于线性化。通过确保$F$的$k$阶Taylor多项式关于高维变量呈线性依赖，扩展ODE系统的线性部分有效地纳入了原始系统的非线性信息。具体而言，(1.1)的数值解按照以下计算框架进行。首先，在每个时间步长$t_n$上，基于数值解点$\mathbf{U}_n$，将一组辅助变量定义为$\mathbf{u}$的最高k次多项式函数。然后推导这些辅助变量的控制方程，得到一个等价的高维ODE系统。通过对此扩展系统应用EIs获得精确解$\mathbf{u}(t_{n+1})$的数值逼近$\mathbf{U}_{n+1}$，随后通过投影步骤将高维解映射回原始状态空间。

我们建立了新类EIs的严格理论分析。一般而言，这些方法对数值解达到$O(h^{k+1})$的误差界，其中$k$为辅助多项式变量的最高次数。此外，对于较大步长$h > c_0\varepsilon$（其中$c_0$为不依赖于$h$和$\varepsilon$的常数），收敛阶提升至$O(\varepsilon h^k)$。将这些结果推广到二阶系统(1.2)时，我们获得了具有改进均匀收敛性的误差估计。一个特别不平凡且关键的认识是：高维系统中按$O(1/\varepsilon)$缩放的线性部分保留了原始矩阵$A$的谱性质。我们通过详细的代数推导严格证明了这一不变性，该结果在建立均匀精度方面起着关键作用。

我们注意到，过去二十年中已开发了众多求解高度振荡微分方程的多尺度方法，以实现均匀精度。例如异质多尺度方法[40]、[41]、多Revolution组合方法[42]、[43]、频闪平均方法[44]、[45]、多尺度时间积分器[12]、[13]、[46]、[47]、嵌套Picard迭代格式[15]、[48]以及双尺度公式格式[22]、[49]。在构造这些数值方法时，许多方法依赖于关于矩阵$A$的结构假设。典型要求包括所有特征值为固定值的整数倍（对应主频率），或者矩阵为对角阵。这些假设有助于解耦系统的振荡模式，从而便于多尺度格式的构造。我们所提方法的一个显著特征是对特征值仅要求零实部的更温和条件。此外，实现高阶均匀精度的数值实现和理论分析仍是一个正在发展的研究领域。我们提出的统一框架——纳入新的线性化技术——为应对这一挑战提供了可行途径。

本文其余部分组织如下。第2节介绍局部线性扩展变量的概念及关联系统，并构造新的EIs。第3节分析高维系统中线性部分的代数性质并建立收敛性分析。第4节给出若干代表性数值实验以验证理论发现。最后一节给出总结性结论。全文中，我们用$C$表示不依赖于$h$和$\varepsilon$的泛常常数。

## 2. 局部线性扩展系统

我们首先将时间$t$视为独立变量，定义增广状态向量$\mathbf{x} = [\mathbf{u}^\top, t]^\top$。原始系统(1.1)改写为：

$$\dot{\mathbf{x}} = \frac{1}{\varepsilon}\mathbf{A}_1\mathbf{x} + \mathbf{f}(\mathbf{x}), \quad \mathbf{x}(0) = \mathbf{x}_{\text{in}},$$

其中

$$\mathbf{A}_1 = (A^{(1)}_{ij}) = \begin{bmatrix} \mathbf{A} & \mathbf{0}_{d \times 1} \\ \mathbf{0}_{1 \times d} & 0 \end{bmatrix}, \quad \mathbf{f}(\mathbf{x}) = \begin{bmatrix} \mathbf{F}(\mathbf{u}, t) \\ 1 \end{bmatrix}, \quad \mathbf{x}_{\text{in}} = \begin{bmatrix} \mathbf{u}_{\text{in}} \\ 0 \end{bmatrix}.$$

此处$\mathbf{0}_{d_1 \times d_2}$表示$d_1 \times d_2$零矩阵。

引入高维变量需要预先定义多元指标和具有关联代数结构的指标集。对于正整数$k$和$d$，定义指标集为：

$$\mathcal{I}[[k]]_{d+1} := \{1, \ldots, d+1\}^{\oplus k},$$

即$d+1$的$k$重包络。

## 3. 关于$A^{[k]}_1$代数性质的初步结果

我们采用以下代数约定。令$\mathcal{P}_j[\mathbf{x}]$表示变量$x_1, \ldots, x_{d+1}$中次数至多为$j$的多项式空间。对于$0 \leq i < j \leq k$，多项式集$\bigcup_{l=i+1}^{j} P[[l, \hat{\mathbf{x}}]]_\mathbf{x}$的线性张成与商空间$\mathcal{P}_j[\mathbf{x}]/\mathcal{P}_i[\mathbf{x}]$自然线性同构。我们将这些同构空间等同，在后续分析中直接记为$\mathcal{P}_j[\mathbf{x}]/\mathcal{P}_i[\mathbf{x}]$。

在实施EIs时，数值格式的性能从根本上由……的性质决定。

## 4. 数值实验

本节通过对(1.2)和(1.1)的数值实验来验证第3节建立的误差界的渐近最优性。采用$k$阶多项式的局部线性扩展指数积分器记为LLEEI$(k+1)$。

**例1** 我们首先验证第2节引入的统一LLEEI框架随着扩展阶数$k$的增加可以实现任意高阶收敛，并表现出关于$\varepsilon$的均匀精度，以及第3节所刻画的逐段行为特征。

## 5. 结论

本文提出了一类基于局部线性扩展技术的EIs。通过采用增维方法，我们对非线性振荡ODE实现了更全面的高阶线性化。具体而言，我们在参考点处引入额外多项式变量来定义局部线性扩展变量，并推导相应的满足系统。对该局部线性扩展系统应用EIs，得到一个统一框架，用于……

---

**CRediT作者贡献声明**

Qi Zhihao：撰写——初稿、可视化、验证、软件、调查、形式分析、概念设计。

Deng Weibing：撰写——审阅与编辑、指导、方法论、形式分析、概念设计。

Zhu Fuhai：形式分析。

**致谢**

我们感谢审稿人富有洞察力和建设性的意见，这些意见显著改进了本文质量。

Qi Zhihao | Deng Weibing | Zhu Fuhai