Oscillatory Differential Equations

Stojimirovic, Tara, and James Bremer. An accelerated frequency-independent solver for oscillatory differential equations. IMA Journal of Numerical Analysis (2026). https://doi.org/10.1093/imanum/draf128

Abstract: …Our method, which operates by constructing a slowly varying phase function representing a basis in the space of solutions of the differential equation, runs in time independent of the frequency of oscillations of the solutions and can be applied to second-order equations whose solutions are oscillatory in some regions and slowly varying in others…

这篇文章[SJ26]考虑求解定义在区间上的二阶常微分方程

要求方程中的函数在时严格正且是光滑的缓慢变化函数(非震荡)，而参数是很大的正数。当时，方程的解是和的线性组合，是以为频率的高频震荡函数。对于一般的系数函数，方程的解通常也是频率与成正比的震荡函数。考虑到一般而言函数的一次振荡需要大概次采样，因此对于上述这种具有超高频率的解函数，常规的数值方法需要非常多的采样点来进行求解，如果使用传统的谱方法，则需要求解一个阶线性系统，使用阶多项式来表示解，而在现实的一些应用中这类方程往往需要成百上千次地求解，例如天文观测中需要的后验Bayes参数估计需要对非常多不同的参数相应的方程进行求解，在类似的这些应用场景下以往的方法并不能满足需求。

下图展示了这类方程的一个例子，该方程对应的解在区域中心快速振荡，在区域两侧变化更平缓。当使用传统方法时，随着频率的增大，为了较好地逼近中心部分的振荡必须使用越来越高阶的多项式构造全局逼近。

振荡微分方程的解，其中。

[SJ26]的核心思想是不直接计算高频振荡的解函数，而是构造原方程的一组解空间：

要求基函数满足原方程，根据这一关系转而求解基函数中的内层函数，这样的称为系统的相函数(phase function)。一旦被计算出，则解可以被表示成

常数与由原方程给定的边界条件确定。

这种方法的主要好处是：尽管解函数非常震荡，但是我们可以找到一个与频率无关的、变化缓慢的光滑相函数来用它导出的基函数表示，得益于这样的的性质足够好，通常只需要远小于的采样点来对它进行插值近似，从而极大程度提高了求解速度。换言之，这一方法的精髓在于通过事先指定解的形式，将一个难以用常规方法表示的极度振荡的函数逼近问题转化为了计算另一个既光滑又非振荡的函数的问题。类似于为了近似某些具有奇异性的函数时，需要构造一组本身含有对应奇异性的正交函数替换经典的正交多项式基进行逼近，我们总是希望根据已有的解函数的结构的信息，将解函数自身的难以处理的结构使用某种转化方法消解掉，这是现代谱方法的一个主要的发展方向，即构造特殊的基函数来让问题得到简化。

基函数的构造

中基函数的构造有明确的动机。让我们从这一最平凡的情形入手，此时方程的解空间由与张成。当不是常数但却只是缓慢变化时，解虽然不像之前那样简单，但却同样具有类似的振荡行为，因此一个自然的想法是寻找形如

形式的解作为解空间的基函数，其中和都是实值函数。求导得

而如果将代回到原方程要求等式成立，就要强制令等式左端的虚部为零，即

两侧同时乘上可得，因此是常数，故

所以满足方程的上述形式的解必须形如。这在物理上对应WKB(Wentzel-Kramers-Brillouin)理论中的作用量守恒或量子力学中的概率流守恒。数学上，WKB方法是求解具有小系数最高阶项的微分方程的最常用的渐近展开方法，常常用来近似近似Schrödinger方程的解。

对于其他类似形式的方程，例如[SJ26]最后提到的

这类方程，其中，由于此时这一点可能对应一个拐点，三角函数不再适合作为解的基础模板，不过我们仍然可以构造类似的广义基函数。以上述方程为例，广义基函数的构造逻辑如下：首先我们选取规范方程，即寻找一个可以反映解的拐点特性的最简单方程，在此处我们考虑的是

接下来，我们假设这一简单方程具有一组线性无关的解与。然后模仿之前的方式，要求原方程的解空间由

张成，这将导出需要满足的广义Kummer方程。求解这一方程就得到了原方程的一组基函数。一个典型的例子是的情形，此时方程的解存在一阶拐点，对应的规范方程自然应当选取Airy方程，按照之前的思路，应该使用Airy函数和来构造基函数，即令

此时通过数值计算就可以避免在原点处的奇异性引发的问题并非常高效地给出原方程的近似解。

算法概述

一旦我们选定了基函数的形式，将基函数代入原方程就可以得到相函数需要满足的方程：

称为Kummer方程。这一非线性方程解空间中的绝大多数解都是高度振荡的，并且对于这一方程我们缺少唯一确定解所需的边界条件，Bremer在之前一篇文章[BR18]中的另一关键发现是：通过区域分割以及Newton迭代，通过欠采样，即在限制区间上采样点的数目，可以强制提取出解空间中唯一的光滑非震荡解。

算法的基本思想如下：首先大致将区间分成一些子区间，这些子区间大致被分为两类，在其中一类区间上原方程的解函数高频振荡，在另一类区间上解函数相对平稳。使用不同的求解器分别计算这两类区间上的。最后以的形式给出基函数并由边界条件确定表出系数。计算的核心在于如何稳定高效地计算两类区间上的。

对于解函数高频振荡的区间，我们考虑计算下述Riccati方程：

它的解与我们需要的之间满足下面的关系：

[SJ26]中证明了当足够大时，使用Chebyshev谱配置法搭配Newton迭代，在选取指定起始函数时无需边界条件可以收敛到唯一的光滑解。这里简要说明为何此处无需边界条件即可进行求解：直接计算可知使用Newton迭代求解Riccati方程时的每次迭代所要计算的线性化方程为

其中，是阶Chebyshev格点上的谱微分矩阵，是在处的Fréchet导数。等价于方程。根据定义

我们选取起始向量为，其中

当满足以下条件：

其中，，则

因此的量级大约为，故当足够大时，矩阵对角占优从而可逆，并且在这种情况下，后续迭代得到的仍然足够大可以保证的可逆性，于是可以收敛到我们想要的唯一的平滑解。事实上，为了求解

我们要求

以保证的谱半径小于1，此时可以使用简单的不动点迭代来完成线性化系统的求解：

另外，条件反过来被用作检验当前区间是否是高频区域的一个判据。

而对于解函数变化较为平缓的区域，我们求解如下Appell方程：

其中与相函数满足以下关系：

作为一个三阶方程，Appell方程需要三个边界条件以保证求解的唯一性，这可以通过在上一步已经得到的相邻的高频区域上的的区间交界处的前几阶导数来指定。之所以需要在较平缓区域使用Appell方程，是因为如果在接近的地方仍然使用上面的方法进行计算会出现数值不稳定的现象，即Riccati方程在拐点附近是病态的，而Appell方程的求解则是在所有的区间上全局稳定的，所以在前述方法不稳定的中低频区域，只需求解Appell方程稳定地得到，再取倒数就可以安全地得到，通过一次谱积分并添加常数使得不同区域交界处连续就可以得到最终的相函数。

简单地来说，这套求解算法采取了“接力”的策略。算法尽可能在所有安全的高频区间内，利用快速且无需边界信息的Newton迭代计算出相位函数的精确值。而当遇到高频方法无法处理的低频或拐点区域时，算法就会提取高频区间的边界值作为初值，切换到更稳定的Appell方程去求解初值或终值问题，从而将相位函数平滑且安全地“延拓”覆盖到整个求解域。