具有递归雅可比自由性质的渐高阶紧致中心格式与柯西-科瓦列夫斯卡亚方法

《Journal of Computational Physics》：Increasinacgly high-order compact central schemes with a recursively Jacobian-free Cauchy-Kovalevskaya method

　　摘要在先前的研究（J. Sci. Comput., 105: 38, 2025）中，我们为加权紧凑型中心差分格式提出了一种无需雅可比矩阵的柯西-科瓦列夫斯卡亚方法。该方法利用泰勒展开根据单元格局部信息计算虚拟点处的通量值，再通过有限差分法求得通量导数，从而避免了繁琐的通量导数张量运算，同时保持了原有加权紧凑型中心差分格式的紧凑性、高精度及稳定性。该算法理论上可实现任意高阶，但在构建更高阶格式时面临两大挑战：其一，高阶导数的有限中心差分公式变得极为复杂；其二，原有算法的CWENO型限幅器在四阶以上格式中效果不佳。本文的主要贡献在于为这类紧凑型中心差分格式提出了一种简化且通用的柯西-科瓦列夫斯卡亚方法，即基于虚拟点的递归有限中心差分公式，用以近似任意高阶的通量导数。尽管这些中心差分仅为二阶，但通过选择合适的虚拟步长，可使数值误差接近机器精度。该方法还保留了原有方法的紧凑性、精度、稳定性及可移植性。为解决第二项挑战，我们用TVB型扰动单元指示器结合斜率限幅器来替代原有的CWENO型限幅器，用于捕捉冲击波，该方法在七阶精度下仍表现良好。大量数值算例证明了该方法的高效性与稳定性。

引言加权紧凑型中心差分格式是一种用于求解双曲守恒定律的高阶数值算法。在每个单元格内，该格式用多项式近似解，并将导数视为独立的自由度。单元格平均值通过基于时空交错网格的经典有限体积格式更新，从而确保相邻控制体之间的界面通量连续，因此无需使用精确的黎曼求解器或近似的黎曼求解器。k阶空间导数是通过单元格顶点处(k?1)阶空间导数的紧凑型中心差分来计算的。为捕捉不连续性，该算法采用了CWENO型限幅器，它有助于保持格式的紧凑性及高精度。所有用于通量计算的时空信息都是通过柯西-科瓦列夫斯卡亚 procedure从空间导数中得到的，这一方法源自柯西与科瓦列夫斯卡亚关于偏微分方程解的存在唯一性定理。在有限差分领域，这一方法也被称为Lax–Wendroff方法。虽然加权紧凑型中心差分格式属于基于交错网格的两步中心差分格式，但其一维中的一阶版本也可重新整理为带数值通量的单步守恒格式。有趣的是，该格式的数值通量恰好是Richtmyer格式与Lax–Friedrichs格式通量的算术平均值，而这正是Toro和Billett所提出的FORCE通量。FORCE格式是Glimm随机选择方法的确定性版本。后续有研究者将一维FORCE格式扩展到多维结构化及非结构化网格上，通过ADER方法实现了空间和时间上的高阶精度。从这个意义上说，加权紧凑型中心差分格式也可视为FORCE格式的高阶扩展版本。

加权紧凑型中心差分格式最突出的优点在于，它能在单个时间步内通过紧凑的格式结构实现均匀的时空高阶精度。实现单时间步内均匀时空高阶精度的关键在于柯西-科瓦列夫斯卡亚procedure，这一方法已被许多其他数值格式广泛采用，如早期的ENO格式、CESE格式、ADER格式、Lax–Wendroff WENO格式以及Lax–Wendroff DG格式等。其中，ADER框架是对这一方法最为成功的应用之一，它通过求解函数值及其空间导数的广义黎曼问题，再利用柯西-科瓦列夫斯卡亚procedure计算时间导数，从而实现了空间和时间上的均匀高阶精度。近年来，ADER格式还被成功与有限体积格式、间断伽辽金有限元格式以及PnPm格式相结合。与流行的多阶龙格-库塔时间离散方法相比，柯西-科瓦列夫斯卡亚procedure每时间步仅需计算一次空间算子，因此在空间算子计算成本较高的情况下具有更好的整体性能。

柯西-科瓦列夫斯卡亚procedure的缺点在于，它不仅依赖于通量函数，还依赖于这些函数对守恒变量的导数。对于一般的非线性系统，高阶导数是高维张量，需要复杂的符号计算，这就导致了该方法的通用性不足。为解决这一问题，Dumbser和Munz提出了自动实现柯西-科瓦列夫斯卡亚procedure的算法。Zhou和Shen则将该技术与加权紧凑型中心差分格式结合，构建了可在统一框架下求解一般守恒定律的任意高阶紧凑型中心差分格式。Montecinos和Balsara在ADER格式的背景下提出了柯西-科瓦列夫斯卡亚procedure的递归公式，简化了代码实现。上述方法均试图简化解析形式的柯西-科瓦列夫斯卡亚procedure的实现，还有些数值方法则试图避免其复杂的符号计算。Demattè等人提出了时间重构方法，无需使用柯西-科瓦列夫斯卡亚procedure即可构建极高阶的ADER格式。Zorío提出了Lax–Wendroff近似泰勒方法，用于改进Qiu和Shu提出的Lax–Wendroff WENO格式。该方法利用高阶中心差分来近似柯西-科瓦列夫斯卡亚方法中的通量导数，简化了实现过程。Carrillo和Parés开发了更为紧凑且稳定性更高的紧凑近似泰勒方法，它通过(2p+1)点格式实现2p阶精度。此外，Carrillo等人还提出了自适应紧凑近似泰勒方法，它能根据局部解的平滑程度调整精度，从而在平滑区域保持高精度，同时在不连续附近抑制振荡。针对平衡方程，人们还开发了平衡型的自适应紧凑近似泰勒方法。Macca等人将高阶紧凑近似泰勒方法与事后多维最优阶数检测方法相结合，形成了CAT-MOOD方法，这类方法在平滑区域具有高精度，在不连续区域几乎无振荡，且数值稳定性优异。Frolkovi?和?áková提出了使用混合时空导数的紧凑隐式格式，这类格式形成了结构较为简单的系统，还可扩展为用于双曲方程的高分辨率TVD格式。Carrillo等人通过将紧凑近似泰勒方法与快速/最优WENO格式结合，降低了其计算成本。Chouchoulis等人为双曲守恒定律设计了无需雅可比矩阵的显式多导数龙格-库塔方法。随后，Chouchoulis和Schütz又为刚性常微分方程设计了无需雅可比矩阵的隐式多导数龙格-库塔方法，通过高阶导数提升了牛顿矩阵的条件数。Moradi等人则为双曲守恒定律设计了无需雅可比矩阵的显式多导数广义线性方法。Lee等人提出了另一种方法，他们在原始皮卡积分公式中引入了无需解析雅可比矩阵和海森矩阵的形式，这种形式类似于在牛顿-克里洛夫迭代中广泛使用的无系统形式，随后他们还通过递归策略将其扩展到四阶，使得实现更加紧凑，效率也更高。最近，有研究者为加权紧凑型中心差分格式提出了无需雅可比矩阵的柯西-科瓦列夫斯卡亚方法，该方法首先通过局部泰勒展开计算虚拟点处的通量，再通过中心有限差分法求得其导数，从而避免了高阶导数张量的复杂符号计算。

在本研究中，我们基于先前的工作，试图构建逐步提高阶数的单阶段紧凑型中心差分格式。然而，该方法在构建更高阶格式时仍面临两大挑战：其一，随着阶数增加，有限中心差分格式需要越来越多的计算点，且各项系数需逐一通过泰勒展开确定；其二，原有算法的CWENO型限幅器在四阶以上格式中效果不佳。为解决第一个问题，我们采用了受Lee等人研究的启发而设计的递归算法，基于虚拟点通过二阶中心差分统一计算任意高阶的通量导数，同时通过选择合适的虚拟步长，使柯西-科瓦列夫斯卡亚procedure的数值误差接近机器精度。将该方法与紧凑型中心差分格式结合，我们构建了可用于求解双曲守恒定律的任意高阶格式，且所有计算都可在同一个程序中完成，我们将该方法称为CCS-RJFCK方法。为解决第二个问题，我们为该方法设计了基于指示器的冲击波捕捉限幅器，用于抑制不连续附近出现的异常振荡。该限幅器实际上是对用于RKDG格式的经典TVB限幅器进行改进而来的，它首先通过改进后的minmod函数识别出所有存在问题的单元格，然后通过TVD型斜率限幅器重建这些单元格中的解。该限幅器与CCS-RJFCK方法的结合，确保了即使在高阶情况下，也能实现清晰且无振荡的不连续面处理。

本文的其余结构如下：第2节概述了基于精确柯西-科瓦列夫斯卡亚procedure的CCS算法框架；第3节详细介绍了基于前述方法的RJFCK方法，并给出了相应的理论分析；第4节阐述了用于冲击波捕捉的基于指示器的限幅器；第5节通过大量数值实验验证了CCS-RJFCK方法的准确性与可靠性；最后，第6节总结了本文的主要研究成果与结论。

任意高阶紧凑型中心差分格式概述关于CCS格式的构建细节可参考[1]和[61]。为便于说明，我们以二维笛卡尔网格为例总结该算法的实现过程。考虑以下二维双曲方程组：?U/?t + ?F(U)/?x + ?G(U)/?y = 0，其中(x,y)属于区域Ω，t属于[0,∞)区间，U表示守恒变量向量，F和G分别代表x方向和y方向的通量函数。

递归无需雅可比矩阵的柯西-科瓦列夫斯卡亚方法用于CCS格式柯西-科瓦列夫斯卡亚procedure使得CCS格式能够在单个时间步内、通过紧凑的格式结构实现均匀的高阶时空精度。然而，导数张量?kF/?Uk和?kG/?Uk的维度为mk+1，因此对于极高阶情况，精确的柯西-科瓦列夫斯卡亚procedure会变得极为复杂，这也就难以将CCS格式推广到更复杂的双曲系统中。

基于指示器的冲击波捕捉限幅器原有的加权紧凑型中心差分格式算法采用了非线性的CWENO型限幅器，该限幅器将CCS格式得到的高阶多项式与若干在特定方向上有偏置的一阶多项式相组合，从而在保持平滑区域的高精度的同时，抑制不连续附近的异常振荡。当设计精度不超过四阶时，该限幅器效果良好，但随着高阶与低阶之间的阶数差距增大，其效果逐渐减弱。

数值结果在本节中，我们通过多个数值算例来验证CCS-RJFCK方法。为便于表述，我们将结合精确柯西-科瓦列夫斯卡亚方法、无需雅可比矩阵的柯西-科瓦列夫斯卡亚方法以及所提方法的k阶CCS格式分别记为CCS-CK-k、CCS-JFCK-k和CCS-RJFCK-k。在所有算例中，无需雅可比矩阵的柯西-科瓦列夫斯卡亚方法的虚拟步长均设置为δx=0.1Δx，δy=0.1Δy，δt=0.1Δt，而RJFCK方法则采用方程(29)中给出的最优步长值。根据表1中的线性稳定性条件，我们进行了相关设置。

结论CCS-RJFCK方法通过在虚拟点基础上递归构建统一的有限中心差分算子，来近似任意高阶的通量导数，从而扩展了基于无需雅可比矩阵的柯西-科瓦列夫斯卡亚方法的框架。因此，该方法允许我们在单个程序中实现任意阶数的CCS-RJFCK格式。尽管有限中心差分算子的形式精度为O(δ2)，但通过合理选择参数，可将计算误差控制在接近机器精度的水平。代码获取：本研究中所生成的定制代码，如需获取，可向相应作者提出合理请求。CRediT作者贡献说明：Run Li：形式分析、方法论、验证、可视化、初稿撰写、概念设计、研究实施。Yuxuan Zhou：形式分析、研究实施、软件开发、方法论、验证。Hua Shen：审稿与编辑、概念设计、形式分析、项目管理、资源协调、监督指导。利益冲突声明：作者声明不存在任何可能影响本文所述工作的已知财务利益或个人关系。Run Li|Yuxuan Zhou|Hua Shen