数学物理学报, 2026, 46(6): 2329-2348

一类线性隐式的保结构松弛方法

蒲浩东, 冉茂华,*

四川师范大学数学科学学院 成都 610066

A Class of Linearly Implicit Structure-Preserving Relaxation Methods

Pu Haodong, Ran Maohua,*

School of Mathematical Sciences, Sichuan Normal University, Chendu 610066

通讯作者: 冉茂华, E-mail: maohuaran@163.com

收稿日期: 2025-08-7   修回日期: 2025-11-25  

基金资助: 四川省科学技术项目(2024NSFSC0441)

Received: 2025-08-7   Revised: 2025-11-25  

Fund supported: Sichuan Science and Technology Program(2024NSFSC0441)

摘要

该文针对一类在工程与物理领域具有重要应用价值的自治常微分方程系统, 提出了一种新型保结构数值算法. 该算法基于一般线性方法的理论框架, 通过引入松弛参数对一般线性方法的系数进行修正, 实现了对系统哈密顿结构、能量守恒律等关键几何特征的精确保持. 所构造的松弛一般线性方法兼具线性隐式和高精度的优势. 理论分析以及针对 Landau-Lifshitz 方程、Kepler 方程、sine-Gordon 方程等物理系统的数值实验表明: 和传统的一般线性方法相比, 松弛修正后的算法在保持原方法绝对稳定域的同时, 计算效率得到了显著提升, 且在结构保持、能量守恒等几何特性方面展现出显著优势. 该方法将松弛技术拓展至一般线性方法框架, 构建了一类适用于自治常微分方程系统的线性隐式保结构算法, 可用于相关物理模型的长时间高精度仿真.

关键词： 常微分方程; 一般线性方法; 松弛技术; 保结构算法

Abstract

This paper proposes a new structure-preserving numerical algorithm for autonomous ordinary differential equation systems with significant applications in engineering and physics. Based on the framework of general linear methods, the algorithm modifies coefficients of conventional methods through relaxation parameters, achieving precise preservation of key geometric features including Hamiltonian structure and energy conservation laws. The constructed relaxed general linear method combines advantages of linear-implicit formulation and high-order accuracy. Theoretical analysis and numerical experiments on typical physical systems-such as the Landau-Lifshitz equation, Kepler equation, and sine-Gordon equation-demonstrate that the relaxed algorithm significantly improves computational efficiency while maintaining the absolute stability region of original methods, and exhibits superior performance in structure preservation and energy conservation. This method extends the relaxation technique to the framework of general linear methods, constructing a class of linearly implicit structure-preserving algorithms suitable for autonomous ordinary differential equation systems, which can be employed for long-time, high-precision simulations of relevant physical models.

Keywords： ordinary differential equations; general linear method; relaxation technique; structure-preserving algorithm

PDF (0KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

蒲浩东, 冉茂华. 一类线性隐式的保结构松弛方法[J]. 数学物理学报, 2026, 46(6): 2329-2348

Pu Haodong, Ran Maohua. A Class of Linearly Implicit Structure-Preserving Relaxation Methods[J]. Acta Mathematica Scientia, 2026, 46(6): 2329-2348

1 引言

考虑自治的常微分系统

$\begin{equation} \begin{cases} y'(t)=f(y(t)),\quad t>0, \\ y(0)=y_0, \end{cases} \end{equation}$

其中 $f:\mathbb{R}^d\to \mathbb{R}^d$ 是一个充分光滑的函数, $y_0\in \mathbb{R}^d$ 是给定初值, 且存在泛函$H:\mathbb{R}^d\to \mathbb{R}$ 使得

$\begin{equation} H(y(t))=H(y(0)),\quad\forall~t\in[T]. \end{equation}$

该自治系统具有广泛的应用场景, 涵盖了诸多知名偏微分方程的空间半离散系统. 如磁动力学领域的经典方程--Landau-Lifshitz 方程, 刻画天体运动的 Kepler 方程以及描述晶体位错运动现象的 sine-Gordon 方程等. 本文的目的是构造一类新的数值算法使得(1.1)的数值解始终处于同一个不变流形上, 即希望数值算法求解(1.1)得到的解序列$\{y_n\}$满足条件

$\begin{equation} y_n\in\mathcal{M}=\{y\big|H(y)=H(y_0)\}. \end{equation}$

这种使得数值解保持微分方程潜在结构的数值算法通常被称为几何数值积分方法或者保结构算法. 这类算法有着显著优势, 它不仅使得数值解能够遵循一定的物理规律, 而且能够有效抑制计算不稳定, 从而获得更可靠的数值解.

在过去几十年里, 保结构算法得到了极大的发展. 众所周知, 传统的 Runge-Kutta (RK) 方法保所有的线性不变量[1], 辛 RK 方法保所有的二次不变量, 但没有 RK 方法能保所有的三次及其以上的多项式不变量[2]. 针对 Hamilton 系统设计保 Hamilton 量的数值格式也被广泛地研究. 如, 冯康等[3]和 Ruth[4] 分别基于生成函数构造了一类辛方法; Marsden 和 West[5] 基于变分原理的离散模拟导出了变分积分子, 发展了离散力学理论; L.Brugnano 等[6]提出了保 Hamilton 量的 Hamilton 边值方法. 上面提到的数值方法有较好的稳定性, 但通常是隐式方法, 且这些方法难以高精度地保持数值解的 Hamilton 量守恒. Simo 等[7]针对一类广泛的机械系统, 构造了能量-动量积分方法. Gonzalez[8] 基于离散方向导数的思想, 引入离散 Hamilton 系统构造了一类有效的守恒算法. 利用这种离散导数的思想, McLachlan 等[9]对更一般的带有 Lyapunov 函数的动力系统的提出了离散梯度方法. 这三种方法虽然都可以精确地保不变量, 但仍然是隐式方法, 难以胜任长时间的数值模拟.

为了克服这一缺陷, Shampine[10] 提出在每一步计算之后, 对获得的数值解加一个小的扰动, 以期使数值解始终处于一个不变流形上. Gear[11] 基于此思想提出了一种保不变量的数值格式, 即投影方法的雏形. 对于 $y_n\in\mathcal{M}$, 投影方法是先用任意一个步长为 $h$ 的单步法 $\Phi_h$ 求解(1.1), 计算 $\widetilde{y}_{n+1}=\Phi_h(y_n)$, 然后令

$\begin{equation} y_{n+1}=\widetilde{y}_{n+1}+\lambda_n H'(\widetilde{y}_{n+1}), \end{equation}$

再通过 $H(y_{n+1})=H(y_0)$ 确定参数 $\lambda_n$. $\lambda_n$ 通常是关于步长 $h$ 的高阶无穷小量, 因此投影方法可以视为通过在 $H'(\widetilde{y}_{n+1})$ 方向上施加微小扰动来促使数值解 $y_{n+1}$ 始终处于同一个不变流形 $\mathcal{M}$. 可以看出如果单步法 $\Phi_h$ 是显式的, 相比于前面提及的隐式保结构算法, 投影方法只需要求解一维的非线性方程而非高维的非线性方程组, 能够显著降低计算成本, 且它与预估阶段所采用的单步法 $\Phi_h$ 具有相同的收敛阶. 基于此思想, Del Buono 等[12]构造了一种保二次不变量的 RK 方法:

针对一个显式的 $s$ 级 RK 方法

$\begin{equation} \begin{cases} Y_{n,1}=y_n,\\ Y_{n,i}=y_n+h\sum^{i-1}_{j=1}\limits a_{ij}f(Y_{n,j}),\quad 2\leq i\leq s, \\ \phi_h(y_n)=y_n+h\sum^s_{j=1}\limits b_jf(Y_{n,j}), \end{cases} \end{equation}$

在增量方向进行扰动可得

$y_{n+1}=y_n+h\lambda_n\sum^s_{j=1}b_jf(Y_{n,j}).$

上式可等价地改写为

$\begin{equation} y_{n+1}=y_n+h\sum^s_{j=1}b_jf(Y_{n,j})+(\lambda_n-1)h\sum^s_{j=1}b_jf(Y_{n,j}).\label{RRK} \end{equation}$

最后通过求解 $H(y_{n+1})=H(y_0)$ 确定参数 $\lambda_n$. 数值结果表明, 参数 $\lambda_n$ 很靠近 $1$. 因此该方法可以视为在增量方向施加微小扰动来促使数值解处于不变流形 $\mathcal{M}$. 这种方法继承了投影方法的优势, 并且具有 RK 方法的仿射不变性质. Ketcheson[13] 称该方法为松弛 RK 方法, 并分析了它的稳定域. 李东方等[14]将松弛技术和隐显 RK 方法结合提出了适用于求解刚性方程的松弛隐显 RK 方法. 谷伟等[15]随后提出了适用于求解高震荡 Hamilton 系统的松弛隐显 RK 方法. Calvo 等[16]对更加一般化的不变量提出了一种嵌入式 RK 方法: 先构造一个比(1.5)低阶的 RK 方法

$\widetilde{\phi}_h(y_n)=y_n+h\sum^s_{j=1}\widetilde{b}_jf(Y_{n,j}),$

然后通过

$\begin{equation} \begin{aligned}\label{ERK} y_{n+1}&=(1-\lambda_n)\phi_h(y_n)+\lambda_n\widetilde{\phi}_h(y_n) \\ &=y_n+h\sum^s_{j=1}b_jf(Y_{n,j})+h\lambda_n\sum^s_{j=1}(\widetilde{b}_j-b_j)f(Y_{n,j}), \end{aligned} \end{equation}$

将低阶的 RK 方法嵌入到(1.5)中; 最后通过 $H(y_{n+1})=H(y_0)$ 确定参数 $\lambda_n$, 以使得数值解处于不变流形上. 可见该方法丰富了扰动方向的选择. Calvo 等[16]还构造了保多个不变量的嵌入式 RK 方法. 容易发现, 无论是嵌入式 RK 方法还是松弛 RK 方法, 都是对 RK 方法某些项的系数施加微小扰动, 且具有类似的表达形式.

基于此, 我们希望将这种对系数施加扰动的思想拓展到更加一般的一般线性方法 (General Linear Methods, GLMs). 与传统保结构方法相比, 本文方法具有以下优势:

· 相较于辛方法和变分积分子: 后两者主要保持辛结构或动量映射, 而本文的方法可直接、精确地保任意指定不变量 (包括非二次和非辛情形);

· 相较于隐式的能量保持方法 (如离散梯度法、Hamilton 边值法): 本文的方法为线性隐式方法, 每步仅需求解一个关于松弛参数的一维非线性方程, 计算成本显著低于求解高维非线性方程组的全隐式方法;

· 相较于投影法: 本文的方法通过将松弛技术用于一般线性方法来实现结构保持. 当 $\lambda_n=0$ 时, 方法即还原为相应的一般线性方法; 通过调节 $\lambda_n$ 则可精确保持不变量. 这种内嵌式的设计赋予了算法更大的普适性与灵活性.

本文其余部分组织如下: 第2节基于一般线性方法, 通过嵌入思想构造了松弛一般线性方法, 并分析了该方法的收敛阶以及松弛参数的存在性; 第3节在特定条件下讨论了所提方法的绝对稳定性, 证明其能保持原基础方法的稳定域特性; 第4节通过数值实验, 验证了算法在收敛精度、能量守恒性能及计算效率方面的有效性; 最后, 第5节对全文工作进行了总结.

2 松弛一般线性方法

本节基于一般线性方法构造保微分方程（1.1）不变量的松弛一般线性方法 (Relaxed General Linear Methods, RGLMs).为了方便, 首先简要回顾一般线性方法的有关基本知识.

2.1 一般线性方法

一般线性方法是多级多值的数值方法, 它涵盖了 RK 方法和线性多步方法.

求解初值问题(1.1)的一般线性方法可描述为:

$\begin{equation} \begin{cases} Y_{i} =\sum_{j=1}^s\limits a_{ij}hF_j+\sum_{j=1}^r\limits u_{ij}y_j^{[n]},\quad i=1,2,\cdots,s,\\ y_{i}^{[n+1]} =\sum_{j=1}^s\limits b_{ij}hF_j+\sum_{j=1}^r\limits v_{ij}y_j^{[n]}, \quad i=1,2,\cdots,r, \end{cases} \end{equation}$

其中 $h$ 是计算步长, $F_i=f(Y_i)$, $y_j^{[n]}$ 是输入值, $y_j^{[n+1]}$ 是输出值. 相应的启动值可以由 RK 方法在内的其他适当方法提供.

一般线性方法(2.1)可以改写为矩阵形式

$ \begin{align*} \left[\begin{array}{c}Y\\y^{[n+1]}\end{array}\right] =\left[\begin{array}{c|c}A\otimes I_d&U\otimes I_d\\\hline B\otimes I_d&V\otimes I_d\end{array}\right] \left[\begin{array}{c}hF\\y^{[n]}\end{array}\right], \end{align*} $

其中 $\otimes$ 为矩阵的 Kronecker 积,

$A=(a_{ij})_{s\times s},~B=(b_{ij})_{r\times s},~U=(u_{ij})_{s\times r},~V=(v_{ij})_{r\times r},$

$I_d\in\mathbb{R}^{d\times d}$ 为单位矩阵,

$ \begin{align*} Y=\left[\begin{array}{c}Y_1\\Y_2\\\vdots\\Y_s\end{array}\right],\quad F=\left[\begin{array}{c}F_1\\F_2\\\vdots\\F_s\end{array}\right],\quad y^{[n]}=\left[\begin{array}{c}y_1^{[n]}\\y_2^{[n]}\\\vdots\\y_r^{[n]}\end{array}\right],\quad y^{[n+1]}=\left[\begin{array}{c}y_1^{[n+1]}\\y_2^{[n+1]}\\\vdots\\y_r^{[n+1]}\end{array}\right]. \end{align*} $

特别地, 当一般线性方法(2.1)以连续 $k$ 步的数值解作为输入和输出时, 即

$\begin{equation} \begin{aligned} y^{[n]}=\big[y_n, y_{n-1},\cdots,y_{n-k+1},hf_n, hf_{n-1},\cdots,hf_{n-k+1}\big]^\top,\\ y^{[n+1]}=\big[y_{n+1}, y_n,\cdots,y_{n-k+2},hf_{n+1}, hf_n,\cdots,hf_{n-k+2}\big]^\top, \end{aligned} \end{equation}$

一般线性方法(2.1)可以表述成如下形式

$\begin{matrix} \begin{cases} &Y_{i} =\sum\limits_{j=1}^sa_{ij}hF_j+\sum\limits_{j=1}^{2k}u_{ij}y_j^{[n]},\quad i=1,2\cdots s,\\ &y_{1}^{[n+1]} =\sum\limits_{j=1}^sb_{1j}hF_j+\sum\limits_{j=1}^{2k}v_{1j}y_j^{[n]},\\ &y_{i}^{[n+1]}=y_{i-1}^{[n]},\quad i=2,\cdots,k,\\ &y_{k+1}^{[n+1]}=hF_s,\\ &y_{i}^{[n+1]}=y_{i-1}^{[n]},\quad i=k+2,\cdots,2k, \end{cases} \end{matrix}$

其中系数满足

$ \begin{align*} a_{sj}=b_{1j},~u_{sj}=v_{1j}. \end{align*} $

特别地, 当 $y^{[n]}=[y_{n},hf_n]$ 时, 一般线性方法(2.3)退化为 RK 方法. 当 $s=1$ 时, 一般线性方法(2.3)将退化为传统的线性 $k$ 步方法. 关于一般线性方法的详细讨论可参见文献[17].

注意到, 基于上面的输入输出设置(2.2), 在一般线性方法(2.3)的输出中, 本质上只需要计算 $y_{n+1}$ (即 $y_{1}^{[n+1]}$), 其余分量均可直接从已知输入 $y^{[n]}$ 中提取. 因此, 简记一般线性方法(2.3)为:

$\begin{matrix} y_{n+1}=y_{1}^{[n+1]}=\sum_{j=1}^sb_{1j}hf(Y_j)+\sum_{j=1}^{2k}v_{1j}y_j^{[n]}. \end{matrix}$

2.2 松弛的一般线性方法

本节将展示基于嵌入思想来构造相应保结构松弛方法的策略.

考虑如下两个一般线性方法:

$\begin{equation} \begin{aligned} \left[\begin{array}{c}\hat{Y}\\\hat{y}^{[n+1]}\end{array}\right] =\left[\begin{array}{c|c}\hat{A}\otimes I_d&\hat{U}\otimes I_d\\\hline\hat{B}\otimes I_d&\hat{V}\otimes I_d\end{array}\right] \left[\begin{array}{c}h\hat{F}\\y^{[n]}\end{array}\right], \end{aligned} \end{equation}$

和

$\begin{equation} \begin{aligned} \left[\begin{array}{c}\tilde{Y}\\\tilde{y}^{[n+1]}\end{array}\right] =\left[\begin{array}{c|c}\tilde{A}\otimes I_d&\tilde{U}\otimes I_d\\\hline \tilde{B}\otimes I_d&\tilde{V}\otimes I_d\end{array}\right] \left[\begin{array}{c}h\tilde{F}\\y^{[n]}\end{array}\right], \end{aligned} \end{equation}$

其中

$ \begin{align*} &y^{[n]}=\big[y_n, y_{n-1},\cdots,y_{n-k+1},hf_n, hf_{n-1},\cdots,hf_{n-k+1}\big]^\top,\\ &\hat{y}^{[n+1]}=\big[\hat{y}_{n+1}, y_n,\cdots,y_{n-k+2},hf_{n+1}, hf_n,\cdots,hf_{n-k+2}\big]^\top,\\ &\tilde{y}^{[n+1]}=\big[\tilde{y}_{n+1}, y_n,\cdots,y_{n-k+2},hf_{n+1}, hf_n,\cdots,hf_{n-k+2}\big]^\top,\\ &\hat{Y}=(\hat{Y}_i)_{s_1\times 1},\hat{F}=\left(f(\hat{Y}_i)\right)_{s_1\times 1},\hat{A}=(\hat{a}_{ij})_{s_1\times s_1},\hat{B}=(\hat{b}_{ij})_{2k\times s_1},\hat{U}=(\hat{u}_{ij})_{s_1\times 2k},\hat{V}=(\hat{v}_{ij})_{2k\times 2k},\\ &\tilde{Y}=(\tilde{Y}_i)_{s_2\times 1},\tilde{F}=\left(f(\tilde{Y}_i)\right)_{s_2\times 1},\tilde{A}=(\tilde{a}_{ij})_{s_2\times s_2},\tilde{B}=(\tilde{b}_{ij})_{2k\times s_2},\tilde{U}=(\tilde{u}_{ij})_{s_2\times 2k},\tilde{V}=(\tilde{v}_{ij})_{2k\times 2k}. \end{align*} $

类似于(2.4), 简记一般线性方法(2.5)和(2.6)分别为

$\begin{matrix}\phi(y^{[n]},h) &= \hat{y}^{[n+1]}_1 = \hat{y}_{n+1} = \sum_{j=1}^{s_1}\hat{b}_{1j}hf(\hat{Y}_j) + \sum_{j=1}^{2k}\hat{v}_{1j}y_j^{[n]}, \end{matrix}$
$\begin{matrix}\psi(y^{[n]},h) &= \tilde{y}^{[n+1]}_1 = \tilde{y}_{n+1} = \sum_{j=1}^{s_2}\tilde{b}_{1j}hf(\tilde{Y}_j) + \sum_{j=1}^{2k}\tilde{v}_{1j}y_j^{[n]}. \label{GLMs2a} \end{matrix}$

联合一般线性方法(2.7)和(2.8), 可构造如下松弛一般线性方法:

$\begin{equation} \begin{aligned} y_{n+1} &=(1-\lambda_n)\phi(y^{[n]},h)+\lambda_n \psi (y^{[n]},h)\\ &=\phi(y^{[n]},h)+\lambda_n\left(\psi (y^{[n]},h)-\phi (y^{[n]},h)\right),\\ \end{aligned} \end{equation}$

其中 $\lambda_n$ 由

$H(y_{n+1})=H(y_0)$

确定.

为了简化描述, 记 $\phi(y^{[n]},h)$ 为 $\phi_n$, $\psi(y^{[n]},h)$ 为 $\psi_n$, 以及

$ \begin{align*} \omega_n: =\omega(y^{[n]},h) =\psi (y^{[n]},h)-\phi (y^{[n]},h). \end{align*} $

则松弛一般线性方法(2.9)可进一步简写为

$\begin{matrix} y_{n+1}=\phi_n+\lambda_n\omega_n. \end{matrix}$

注 2.1 从形式上看, 在松弛一般线性方法(2.9)的构造中, 存在方法(2.7)和方法(2.8)有相同长度的输入和输出的限制. 但实际上该约束不是本质的, 因为两个方法中长的输入和输出总能包括短的输入和输出, 不足部分视方法的相应系数为 $0$ 即可.

下面讨论 $\lambda_n$ 的存在性. 在这之前, 我们给出如下引理.

引理 2.1 设 $a(x)$ 是 $\mathbb{R}$ 到 $\mathbb{R}$ 的连续函数,

$f(x)=a(x)x^2+bx+c,$

其中 $b,c$ 是常数且 $b\neq 0$. 若存在常数 $M$ 对任意 $x\in\mathbb{R}$ 满足

$\begin{matrix} M>b^2-4a(x)c\geqslant 0, \end{matrix}$

则 $f(x)=0$ 有解.

证 当 $c=0$ 时, $f(0)=c=0$, 此时 $f(x)=0$ 有零解.

当 $c\neq 0$ 时, 有

$\frac{4c}{x^2}f(x)=4c^2\left(\frac{1}{x}+\frac{b}{2c}\right)^2-\left(b^2-4a(x)c\right).$

令

$g(x)=4c^2\left(\frac{1}{x}+\frac{b}{2c}\right)^2-\left(b^2-4a(x)c\right),$

则 $f(x)=0$ 有解等价于 $g(x)=0$ 有解. 注意到假设(2.11), 有

$\lim\limits_{x\to 0} g(x)= +\infty,\quad g(-\frac{2c}{b}) \leqslant 0.$

由介值定理可知 $f(x)=0$ 有解. 证毕.

接下来讨论 $\lambda_n$ 的存在性和松弛一般线性方法的收敛阶.为了方便叙述, 记 $\frac{\omega_n}{\left\|\omega_n\right\| }=d_n$.

定理 2.1 假设方法(2.5)是 $p$ 阶的, $H(y)$ 二阶可微且满足 $\left|H''(y)\right|<M$, $M$ 是常数. 记 $ H'(\phi_n)^\top d_n/h^r=B(y^{[n]},h)$, 其中整数 $r\ge0$ 且 $2r<p+1$. 如果当 $h\to 0$ 时, $B(y^{[n]},h)\to B(y^{[n]},0)$ 且 $B(y^{[n]},0)\neq 0$. 那么存在 $h^*>0$, 使得对任意的 $h\in[0,h^*]$, 存在 $\lambda_n$ 使得 $H(\phi_{n}+\lambda_n\omega_n)=H(y_0)$, 且相应的松弛一般线性方法(2.10)是 $p-r$ 阶的.

特别地当 $\omega_n=H'(\phi_n)$ 时, 方法是 $p$ 阶的; 当 $\omega_n=\phi_n-y_n$ 时, 方法是 $p-1$ 阶的.

证 记 $\lambda_n\left\|\omega_n\right\|=\mu$, 令

$g(\mu)=H\left(\phi_{n}+\lambda_n\omega_n\right)-H(y_0)=H\left(\phi_{n}+\mu d_n\right)-H(y_0).$

将 $g(\mu)$ 在 $\mu=0$ 处进行Taylor展开可得

$g(\mu)=H(\phi_{n})-H(y_0)+\mu H'(\phi_n)^\top d_n+\frac{\mu^2}{2}d^\top_nH''(\phi_n+\xi d_n)d_n,$

其中 $\xi\in (0,\mu)$ 是关于 $\mu$ 的函数. 由中值定理知

$H(\phi_{n})-H\left(y(0)\right) =H(\phi_n)-H\left(y(t_n)\right) = H'\left(\phi_{n}+\xi_1\left(y(t_n)-\phi_{n}\right)\right)^\top\left(\phi_n-y(t_n)\right)=\mathcal{O}(h^{p+1}),$

其中 $\xi_1\in(0,1)$.

记 $H(\phi_{n})-H\left(y(0)\right)=a_nh^{p+1}$, $a_n$ 是一个有界量. 由于 $y(0)=y_0$, 则有

$g(\mu)=a_nh^{p+1}+h^{r}\mu B(y^{[n]},h)+\frac{\mu^2}{2}d_n^\top H''(\phi_n+\xi d_n)d_n.$

根据引理2.1, 函数 $g(\mu)$ 是否有解只需判断其系数是否满足引理2.1的条件. 因为 $B(y^{[n]},0)\neq 0$, 所以当 $h$ 充分小时, $\left(B(y^{[n]},h)\right)^2>0$ 且有界. 又因为 $2r<p+1$, 所以当 $h$ 充分小时, 有

$2\left(h^rB(y^{[n]},0)\right)^2>h^{2r}\left(B(y^{[n]},h)\right)^2-2h^{p+1}a_nd^\top_nH''(\phi_n+\xi d_n)d_n\ge 0.$

因此, 由引理2.1知 $g(\mu)=0$ 有解, 即存在满足条件的 $\lambda_n$.

当 $d^\top_nH''(\phi_n+\xi d_n)d_n=0$ 时, $g(\mu)=0$ 的解为

$\mu=-h^{p-r+1}a_n/B(y^{[n]},h)=\mathcal{O}(h^{p-r+1}).$

当 $d^\top_nH''(\phi_n+\xi d_n)d_n\neq 0$ 时, 若 $B(y^{[n]},h)>0$, 由求根公式得

$\mu=\frac{2a_nh^{p-r+1}}{-B(y^{[n]},h)-\sqrt{\left(B(y^{[n]},h)\right)^2-2h^{p-2r+1}a_nd^\top_nH''(\phi_n+\xi d_n)d_n}}=\mathcal{O}(h^{p-r+1})$

是 $g(\mu)=0$ 的解; 若 $B(y^{[n]},h)<0$, 则

$\mu=\frac{2a_nh^{p-r+1}}{-B(y^{[n]},h)+\sqrt{\left(B(y^{[n]},h)\right)^2-2h^{p-2r+1}a_nd^\top_nH''(\phi_n+\xi d_n)d_n}}=\mathcal{O}(h^{p-r+1})$

是 $g(\mu)=0$ 的解. 故 $\lambda_n\omega_n=\mu d_n=\mathcal{O}(h^{p-r+1})$. 进一步, 由松弛一般线性方法(2.10)的定义可知, $y_{n+1}$ 等价于 $\phi_n$ 加上一个 $p-r+1$ 阶的扰动项. 这意味着, 相应的松弛一般线性方法(2.10)是 $p-r$ 阶的.

特别地, 对于正交投影 $\omega_n=H'(\phi_n)$, 有 $ H'(\phi_n)^\top d_n=\left\| H'(\phi_n)\right\|$. 这意味着 $r=0$, 即正交投影不会降低精度的阶.

当 $\omega_n=\phi_n(y^{[n]},h)-y_n$ 时, 有

$\omega_n=\phi_n(y^{[n]},h)-y_n=y(t_{n+1})-y(t_n)+\mathcal{O}(h^{p+1})=hy'(\xi)+\mathcal{O}(h^{p+1}),$

其中 $\xi\in(t_n,t_{n+1})$. 于是有

$ \begin{align*} H'(\phi_n)^\top d_n &=\left(H'(\phi_n)- H'\left(y(t_{n+1})\right)\right)^\top d_n\\ & + H'\left(y(t_{n+1})\right)^\top\left(d_n-\frac{y'(t_{n+1})}{\left\|y'(\xi)\right\|}\right)+ H'\left(y(t_{n+1})\right)^\top\frac{y'(t_{n+1})}{\left\|y'(\xi)\right\|}. \end{align*} $

因为 $H'(\phi_n)- H'\left(y(t_{n+1})\right)=\mathcal{O}(h^{p+1}),~H'\left(y(t_{n+1})\right)^\top y'(t_{n+1})=0$, 所以

$H'(\phi_n)^\top d_n = H'\left(y(t_{n+1})\right)^\top\left(d_n-\frac{y'(t_{n+1})}{\left\|y'(\xi)\right\|}\right)+\mathcal{O}(h^{p+1}).$

由于

$ \begin{align*} d_n&=\frac{\omega_n}{\left\|\omega_n\right\|}=\frac{y'(\xi)}{\left\|y'(\xi)+\mathcal{O}(h^p)\right\|}+\mathcal{O}(h^p)\\ &=\frac{y'(\xi)}{\left\|y'(\xi)\right\|}+y'(\xi)\frac{\left\|y'(\xi)\right\|-\left\|y'(\xi)+\mathcal{O}(h^p)\right\|}{\left\|y'(\xi)\right\|\left\|y'(\xi)+\mathcal{O}(h^p)\right\|}+\mathcal{O}(h^p)\\ &=\frac{y'(\xi)}{\left\|y'(\xi)\right\|}+\mathcal{O}(h^p), \end{align*} $

把 $d_n$ 代入 $H'(\phi_n)^\top d_n$ 得到

$H'(\phi_n)^\top d_n = H'(y(t_{n+1}))^\top\frac{y'(\xi)-y'(t_{n+1})}{\left\|y'(\xi)\right\|}+\mathcal{O}(h^p).$

因为

$ \begin{align*} y'(\xi)-y'(t_{n+1})&=\frac{y(t_{n+1})-y(t_n)}{h}-y'(t_{n+1})\\ &=\frac{hy'(t_{n+1})-h^2y''(\xi_1)/2}{h}-y'(t_{n+1})\\ &=-\frac{hy''(\xi_1)}{2}, \end{align*} $

所以 $ H'(\phi_n)^\top d_n=-h H'(y(t_{n+1}))^\top \frac{y''(\xi_1)}{2\left\|y'(\xi)\right\|}+\mathcal{O}(h^p)$. 这意味着 $r=1$, 所以收敛阶会降低一阶. 证毕.

注 2.2 当基础方法 $\phi_n$ 为 RK 方法, 且取 $\omega_n=\phi_n-y_n$ 时, 该松弛一般线性方法将退化为文献[13] 中的松弛RK方法.

上述定理给出了 $\omega_n$ 满足什么条件时, 方法能够保持不变量. 如果方法(2.7)是 $p$ 阶方法, 方法(2.8)是 $q$ 阶方法, 不妨假设 $p>q$, 记

$y^{[n]}_*=\big[y(t_n), y(t_{n-1}),\cdots y(t_{n-k+1}),hf(t_n), hf(t_{n-1}),\cdots hf(t_{n-k+1})\big]^\top,$

则 $\omega_n=\phi(y^{[n]}_*,h)-\psi(y^{[n]}_*,h)$ 是关于 $h$ 的 $q$ 阶无穷小量.

为了便于后面稳定性的分析, 我们给出如下定理.

定理 2.2 在定理2.1假设下, 如果 $\phi_n(y^{[n]},h)$ 和 $\phi_n(y^{[n]},h)+\lambda_n\omega_n(y^{[n]},h)$ 都是 $p$ 阶方法, $\omega_n(y^{[n]}_*,h)$ 是关于 $h$ 的 $q$ 阶无穷小量且 $q<p$, 那么 $\lambda_n=\mathcal{O}(h^{p-q})$.

证 因为 $\phi_n(y^{[n]},h)+\lambda_n\omega_n(y^{[n]},h)$ 是 $p$ 阶方法, 所以

$ \begin{align*} \phi_n(y^{[n]},h)+\lambda_n\omega_n(y^{[n]},h)-\left(\phi_n(y^{[n]}_*,h)+\lambda_n\omega_n(y^{[n]}_*,h)\right) = \mathcal{O}(h^{p+1}). \end{align*} $

即 $ \begin{align*} \phi_n(y^{[n]},h)-\phi_n(y^{[n]}_*,h)+\lambda_n\left(\omega_n(y^{[n]},h)-\omega_n(y^{[n]}_*,h)\right)=\mathcal{O}(h^{p+1}). \end{align*} $

又因为 $\phi_n(y^{[n]},h)$ 也是 $p$ 阶方法, 进而有

$ \begin{align*} \mathcal{O}(h^{p+1})+\lambda_n(\omega_n(y^{[n]},h)-\omega_n(y^{[n]}_*,h))=\mathcal{O}(h^{p+1}). \end{align*} $

注意到 $\omega_n(y^{[n]}_*,h)$ 是关于 $h$ 的 $q$ 阶无穷小量且 $q<p$, 进一步有

$ \begin{align*} \mathcal{O}(h^{p+1})+\lambda_n\mathcal{O}(h^{q+1})=\mathcal{O}(h^{p+1}). \end{align*} $

即 $\lambda_n=\mathcal{O}(h^{p-q})$. 证毕.

3 稳定性分析

直接分析松弛一般线性方法(2.9)的绝对稳定性是困难的, 本节考虑一类满足如下条件3.1的松弛一般线性方法的绝对稳定性.

条件 3.1 假设构造松弛一般线性方法(2.9)的两个基础线性方法(2.5)和(2.6)是显式的, 且其系数满足 $s_1=s_2,~\hat{a}_{ij}=\tilde{a}_{ij},~~i=1,\cdots,s_1-1,~j=1,\cdots,2k.$

即矩阵 $\hat{A}$ 和矩阵 $\tilde{A}$ 的前 $s_1-1$ 行元素相同.

注 3.1 条件3.1并不严苛. 如, 对于给定的显式一般线性方法(2.5), 满足假设条件3.1的一般线性方法(2.6)具有如下形式:

$ \begin{align*} \begin{cases} &\tilde{Y}_i=\sum\limits_{j=1}^{s_1}\hat{a}_{ij}hF_j+\sum\limits_{j=1}^{2k}\hat{u}_{ij}y_j^{[n]},~i=1,2,\cdots s_1\!-\!1,\\ &\tilde{Y}_{s_1}=\sum_{j=1}^{2k}\limits \tilde{v}_{s_1j}y_j^{[n]},\\ &\tilde{y}_1^{[n+1]}=\sum_{j=1}^{2k}\limits \tilde{v}_{s_1j}y_j^{[n]},\\ &\tilde{y}_{i}^{[n+1]}=y_{i-1}^{[n]},\quad i=2,\cdots,k,\\ &\tilde{y}_{k+1}^{[n+1]}=hf(\tilde{Y}_{s_1}),\\ &\tilde{y}_{i}^{[n+1]}=y_{i-1}^{[n]},\quad i=k+2,\cdots,2k. \end{cases} \end{align*} $

事实上, 显式的线性 $k$ 步方法均可写成上述形式. 因为, 针对显式的线性 $k$ 步方法而言, 只有 $\tilde{Y}_{s_1}$ 才参与 $\tilde{y}^{[n+1]}$ 的计算, 其他 $\tilde{Y}_{i}(i=1,2,\cdots,s_1\!-\!1)$ 并不参与 $\tilde{y}^{[n+1]}$ 的计算. 因此可以任意选取计算 $\tilde{Y}_{i}(i=1,2,\cdots,s_1\!-\!1)$ 的方法系数, 自然也包括和方法(2.5)一样的系数.

满足假设条件3.1的松弛一般线性方法可以表述为:

$\begin{matrix} \left[\begin{array}{c}Y\\y^{[n+1]}\end{array}\right]=\left[\begin{array}{c|c}(\hat{A}+\lambda_n A')\otimes I_d&(\hat{U}+\lambda_n U')\otimes I_d\\\hline (\hat{B}+\lambda_n B')\otimes I_d&(\hat{V}+\lambda_n V')\otimes I_d\end{array}\right]\left[\begin{array}{c}hF\\y^{[n]}\end{array}\right], \end{matrix}$

其中

$A'=\tilde{A}-\hat{A},~B'=\tilde{B}-\hat{B},~U'=\tilde{U}-\hat{U},~V'=\tilde{V}-\hat{V}.$

考虑测试方程

$\begin{matrix} y'(t)=qy(t), {\rm Re}(q)<0. \end{matrix}$

记 $z=qh$, 应用松弛一般线性方法(3.1)到测试方程(3.2)得

$\begin{matrix} \left[\begin{array}{c}Y\\y^{[n+1]}\end{array}\right]=\left[\begin{array}{c|c}\hat{A}+\lambda_n A'&\hat{U}+\lambda_n U'\\\hline\hat{B}+\lambda_n B'&\hat{V}+\lambda_n V'\end{array}\right]\left[\begin{array}{c}zY\\y^{[n]}\end{array}\right]. \end{matrix}$

整理可得

$\begin{matrix} y^{[n+1]}=\left(\hat{V}+\lambda_n V' +z(\hat{B}+\lambda_nB')\left(I-z(\hat{A}+\lambda_nA')\right)^{-1}(\hat{U}+\lambda_nU')\right)y^{[n]}. \end{matrix}$

称

$\begin{matrix} M_{\lambda_n}(z)=\hat{V}+\lambda_n V' +z(\hat{B}+\lambda_nB')\left(I-z(\hat{A}+\lambda_nA')\right)^{-1}(\hat{U}+\lambda_nU') \end{matrix}$

为松弛一般线性方法(3.1)的稳定矩阵. 特别地, 当 $\lambda_n\equiv0$ 时, 松弛一般线性方法(3.1)将退化为一般线性方法(2.5), 相应的稳定矩阵也将退化为

$\begin{matrix} M_0(z)=\hat{V}+z\hat{B}(I-z\hat{A})^{-1}\hat{U}. \end{matrix}$

借助于稳定矩阵, 类似于文献[17]中对一般线性方法的绝对稳定域的定义, 给出如下关于松弛一般线性方法的绝对稳定域定义.

定义 3.1 对于松弛一般线性方法(3.1), 定义其绝对稳定域为

$\begin{matrix} \mathbb{D}=\left\{z\in\mathbb{C}\big|\sup^{\infty}_{n=1}\left\|\prod^{n}_{i=1}M_{\lambda_i}(z)\right\|<\infty\right\}. \end{matrix}$

不难验证, 当 $\lambda_n\equiv 0$ 时, 该定义将随之退化为文献[17]中关于一般线性方法的绝对稳定域定义.

注意到条件3.1, 矩阵 $A'$ 的前 $s-1$ 行元素全为零, 即 $A'$ 的秩为 $1$. 这意味着 $A'$ 可表示为 $A'=uv'$, 其中 $u,v$ 是 $s$ 维列向量. 因此, 有

$\begin{matrix} \left(I-z(\hat{A}+\lambda_nA')\right)^{-1}=(I-z\hat{A}-z\lambda_nuv')^{-1}=(I-z\hat{A})^{-1}+\alpha\lambda_nA_1,\label{eq:3.6} \end{matrix}$

其中

$\begin{matrix} \alpha=\frac{z}{1-z\lambda_nv'(I-z\hat{A})^{-1}u},~~ A_1(z)=(I-z\hat{A})^{-1}uv'(I-z\hat{A})^{-1}. \end{matrix}$

根据(3.5), 进一步有

$\begin{matrix} M_{\lambda_n}(z)&=\hat{V}+\lambda_nV'+z(\hat{B}+\lambda_nB')\left((I-z\hat{A})^{-1}+\lambda_n\alpha A_1\right)(\hat{U}+\lambda_nU')\nonumber\\ &=M_0(z)+\lambda_nN_1(\alpha,z)+\lambda_n^2N_2(\alpha,z)+\lambda_n^3N_3(\alpha,z),\label{eq:3.8} \end{matrix}$

其中

$ \begin{align*} &N_1(\alpha,z)=zB'(I-z\hat{A})^{-1}\hat{U}+z\hat{B}(I-z\hat{A})^{-1}U'+\alpha z\hat{B}A_1\hat{U}\triangleq N_{11}(z)+\alpha N_{12}(z),\\ &N_2(\alpha,z)=zB'(I-z\hat{A})^{-1}U'+\alpha(zB'A_1\hat{U}+z\hat{B}A_1U')\triangleq N_{21}(z)+\alpha N_{22}(z),\\ &N_3(\alpha,z)=\alpha zB'A_1U'\triangleq N_{31}(z)+\alpha N_{32}(z). \end{align*} $

从(3.10)可以看出, 松弛一般线性方法的稳定矩阵是对应一般线性方法稳定矩阵的扰动. 下面分析松弛一般线性方法(3.1)的绝对稳定性.

定理 3.1 若对应一般线性方法(2.5)的稳定矩阵满足 $\left\|M_0(z)\right\|<1$, 则相应的松弛一般线性方法(3.1)也是绝对稳定的.

证 注意到

$\left\|M_{\lambda_n}(z)\right\|\le \left\|M_0(z)\right\|+\sum^3_{i=1}|\lambda_n|^i\left(\left\| N_{i1}(z)\right\|+|\alpha|\left\| N_{i2}(z)\right\|\right).$

由(3.9)可知, 当 $|\lambda_n|<\frac{1}{2|zv'(I-zA)^{-1}u|}$ 时, 有 $|\alpha|<2|z|$. 此时, 有

$\left\|M_{\lambda_n}(z)\right\|\le \left\|M_0(z)\right\|+\sum^3_{i=1}|\lambda_n|^i\left(\left\| N_{i1}(z)\right\|+2\left\|zN_{i2}(z)\right\|\right).$

因为当 $\lambda_n=0$ 时, 有

$\left\|M_{\lambda_n}(z)\right\|=\left\|M_0(z)\right\|<1.$

由定理2.2知 $\lambda_n$ 是关于步长 $h$ 的高阶无穷小量. 所以当 $h$ 充分小时, 有

$\sum^3_{i=1}|\lambda_n|^i\left(\left\| N_{i1}\right\|+2\left\|zN_{i2}\right\|\right)<\frac{1-\left\|M_0(z)\right\|}{2}.$

所以有

$\left\|M_{\lambda_n}(z)\right\|<\frac{1+\left\|M_0(z)\right\|}{2}<1.$

即

$\left\|\prod^{n}_{i=1}M_{\lambda_i}(z)\right\|<1.$

所以相应的松弛一般线性方法(3.1)绝对稳定.

证毕.

图1给出了后续数值算例涉及的几个松弛一般线性方法在不同的松弛参数 $\lambda$ 下的绝对稳定域. 容易观察到, 当 $|\lambda|\to 0$ 时, 松弛一般线性方法和对应的基本方法 (i.e., $\lambda\!=\!0$) 有相同的绝对稳定域.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1   不同松弛方法的绝对稳定域


4 数值实验

本节将采用若干代表性模型对所提数值算法的收敛率及守恒性能进行验证. 为便于比较, 首先给出计算过程中所涉及的一些相关指标的定义.

·误差 : 定义在时间区间 $[T]$ 上的全局误差:

$ \begin{align*} E(h) = \max_{0 \leq n \leq N} \| y_n - y(t_n) \|_\infty, \end{align*} $

其中 $y(t_n)$ 为精确解, ${y_n}$ 为数值解, $N = T/h$ 为总步数.

·收敛阶 : 定义收敛阶:

$ \begin{align*} p=\log_2 \left( \frac{E(h)}{E(h/2)} \right). \end{align*} $

·守恒误差 : 定义不变量 $H$ 的相对误差为:

$ \begin{align*} \frac{\|H(y_n) - H(y_0)\|}{\|H(y_0)\|}, \end{align*} $

用以评估算法的保结构性能.

例 4.1 考虑磁动力学领域的经典方程-Landau-Lifshitz 方程[16]:

$\begin{matrix} \begin{cases} y'=H_{\mathrm{eff}}\times y+\lambda y\times (H_{\mathrm{eff}}\times y),\\ y(0)=(\sqrt{6}/4,-\sqrt{6}/4,1/2)^\top, \end{cases} \end{matrix}$

其中 $\lambda=1/20.1,~H_{\mathrm{eff}}=(1,0,0)^\top$, 对应的精确解为:

$y(t)=\left(\frac{a(t)}{b(t)},\frac{2}{b(t)}(-\frac{\sqrt{6}}{4}\cos t-\frac12\sin t),\frac{2}{b(t)}(-\frac{\sqrt{6}}{4}\sin t+\frac12\cos t)\right)^\top,$

其中

$ \begin{align*} a(t)={\rm e}^{\lambda t}(1+\frac{\sqrt{6}}{4})-{\rm e}^{-\lambda t}(1-\frac{\sqrt{6}}{4}), \quad b(t)={\rm e}^{\lambda t}(1+\frac{\sqrt{6}}{4})+{\rm e}^{-\lambda t}(1-\frac{\sqrt{6}}{4}). \end{align*} $

该方程主要用于描述磁性材料中磁化矢量在有效磁场作用下的动态行为.

注意到 $y'(t)$ 与 $y(t)$ 正交, 因此 $H(y)=y^\top y$ 是二次不变量. 显然, 有 $H(y(0))=1$, 这意味着方程(4.1)的解位于单位球面上.

为了便于比较, 我们考虑文献[18]中的五阶混合多步方法 (HMM):

$\begin{matrix} \begin{cases} &Y_1=-\frac{529}{3375}y^{[n-1]}_1+\frac{3904}{3375}y^{[n-1]}_2+\frac{4232}{3375}y^{[n-1]}_3+\frac{1472}{3375}y^{[n-1]}_4,\\ &Y_2=\frac{189}{92}hf(Y_1)+\frac{152}{25}y^{[n-1]}_1-\frac{127}{25}y^{[n-1]}_2+-\frac{419}{100}y^{[n-1]}_3-\frac{1118}{575}y^{[n-1]}_4,\\ &Y_3=\frac{3375}{5152}hf(Y_1)+\frac{25}{168}hf(Y_2)+y^{[n-1]}_1+\frac{19}{96}y^{[n-1]}_3-\frac{1}{552}y^{[n-1]}_4,\\ &y^{[n]}_1=\frac{3375}{5152}hf(Y_1)+\frac{25}{168}hf(Y_2)+y^{[n-1]}_1+\frac{19}{96}y^{[n-1]}_3-\frac{1}{552}y^{[n-1]}_4,\\ &y^{[n]}_2=y^{[n-1]}_1,\\ &y^{[n]}_3=hf(Y_3),\\ &y^{[n]}_4=y^{[n-1]}_3, \end{cases} \end{matrix}$

其中

$ \begin{align*} y^{[n-1]}=\big[y_{n-1},y_{n-2},hf_{n-1},hf_{n-2}\big]^\top,~~ y^{[n]}=\big[y_n,y_{n-1},hf_n,hf_{n-1}\big]^\top. \end{align*} $

将向前 Euler 方法:

$\begin{matrix} \tilde{y}^{[n]}_1=y^{[n-1]}_1+y^{[n-1]}_3 \end{matrix}$

嵌入到该混合多步方法中得到松弛混合多步方法 (RHMM):

$\begin{matrix} y_n=y^{[n]}_1+\lambda_n(\tilde{y}^{[n]}_1-y^{[n]}_1). \end{matrix}$

· 注意到不变量 $H(y) = y^\top y$, 因此其梯度 $H'(y) = 2y$, Hessian 矩阵 $H''(y) = 2I$. 显然, 有 $\|H''(y)\| = 2$. 这意味着满足定理2.1的光滑性假设和有界性假设.

· 由于该算例的精确解满足 $y(t)^\top y(t) = 1$, 因此有 $y^\top y' = 0$ 及 $y^\top y'' = -\|y'\|^2$. 记方法(4.2)为 $\phi_n$, 向前 Euler 方法为 $\psi_n$. 由 Taylor 展开可得

$\begin{matrix} H'(\phi_n)^\top d_n = \frac{-2 y(t_n)^\top y''(t_n)}{\|y''(t_n)\|}+ O(h) =\frac{2 \|y'(t_n)\|^2}{\|y''(t_n)\|} + O(h). \end{matrix}$

注意到方程(4.1) 的向量场在 $y$ 平行于 $H_{\mathrm{eff}}$ 时为零 (即该方向为孤立平衡点), 而初值 $y(0)$ 显然不平行于常向量 $H_{\mathrm{eff}} = (1,0,0)^\top$. 由解的唯一性可知, $y(t)$ 始终不与 $H_{\mathrm{eff}}$ 平行, 故 $H_{\mathrm{eff}} \times y(t) \neq 0$. 此外, 由方程(4.1)知

$\|y'(t)\|^2 = \|H_{\mathrm{eff}} \times y(t)\|^2 + \lambda^2 \|y(t) \times (H_{\mathrm{eff}} \times y(t))\|^2> 0.$

由 $\|y'(t)\|\neq0$ 及 $y^\top y'' = -\|y'\|^2$ 可知 $\|y''\|\neq 0$, 因此当 $h \to 0$ 时, 有

$\begin{matrix} B(y^{[n]},0) = 2\|y'(t_n)\|^2/\|y''(t_n)\| \neq 0. \end{matrix}$

这意味着, 定理2.1的另一个假设也成立.

根据定理2.1, $r$ 由

$ \begin{align*} \frac{H'(\phi_n)^\top d_n}{h^r} \to B(y^{[n]},0) \ne 0 \quad (h \to 0), \end{align*} $

确定. 注意到(4.5)-(4.6), 可得 $r = 0$. 又因为 HMM 方法的阶数 $p = 5$, 根据定理2.1, 松弛方法的收敛阶不低于 $5$ 阶.

分别用 RHMM 和 HMM 方法在 $[0,10^4]$ 的时间区间内求解方程(4.1), 相应的计算结果详见表1. 容易观察到, RHMM 方法的数值收敛阶位于 $4.9$ 至 $7.2$ 之间. 即使针对大步长 $h=5/3$, 对应的收敛阶也接近 $5$. 针对更小的步长, 甚至有收敛阶 $>6$ 的超收敛现象.

表1   HMM 方法和 RHMM 方法对应的误差和收敛阶

新窗口打开| 下载CSV


为评估 RHMM 方法的保结构性能, 图2展示了 $h=5/6$ 时不变量 $H$ 的相对误差. 可见 RHMM 方法能精确保持不变量 (数值解满足 $y_n^\top y_n = 1$), 而 HMM 无法保持该二次不变量. 这些数值结果支持了理论分析结果的正确性.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2   $h=5/6$ 时 Landau-Lifshitz 方程不变量的相对误差演化图


例 4.2 考虑描述两个相互吸引的物体运动过程的 Kepler 方程. 设第一个物体固定在原点, 记第二个物体的位置和速度分别为 $q=(q_1, q_2)^\top,p=(p_1, p_2)^\top$, 则 $p,q$ 满足 Hamilton 系统

$\begin{matrix} \left[\begin{array}{c}p'\\q'\end{array}\right]=J\nabla H,\quad J=\left[\begin{array}{cc}0&-I\\I&0\end{array}\right], \end{matrix}$

其中 Hamilton 量

$H(p,q)=\frac12\left\|p\right\|^2_2-\frac{1}{\left\|q\right\|_2}.$

参照文献[19], 取初值条件

$ \begin{align*} q(0)=(1-e,0)^\top\quad p(0)=(0,\sqrt{\frac{1+e}{1-e}})^\top. \end{align*} $

当离心率 $e\in(0,1)$ 时, 第二个物体的运动轨迹是以第一个物体的位置为焦点的椭圆. $q$ 的精确解为

$ \begin{align*} q_1(t)=\cos (E)-e,\quad q_2(t)=\sqrt{1-e^2}\sin{E}, \end{align*} $

其中参数 $E$ 由 Kepler 公式 $t=E-e\sin(E)$ 决定.

类似地, 考虑文献[20]中给出的三阶两步 RK 方法 (TSRK):

$\begin{matrix} \begin{cases} &Y_1=y^{[n]}_1+\frac 12y^{[n]}_3,\\ &Y_2=y^{[n]}_2+\frac 12y^{[n]}_4,\\ &Y_3=\frac56 hf(Y_1)-\frac56 hf(Y_2)+y^{[n]}_1+\frac23y^{[n]}_3+\frac13y^{[n]}_4,\\ &y^{[n+1]}_1=\frac56 hf(Y_1)-\frac56 hf(Y_2)+y^{[n]}_1+\frac23y^{[n]}_3+\frac13y^{[n]}_4,\\ &y^{[n+1]}_2=y^{[n]}_1,\\ &y^{[n+1]}_3=hf(Y_3),\\ &y^{[n+1]}_4=y^{[n]}_3, \end{cases} \end{matrix}$

其中

$ \begin{align*} y^{[n]}=\big[y_{n},y_{n-1},hf_{n},hf_{n-1}\big]^\top,~~ y^{[n+1]}=\big[y_{n+1},y_{n},hf_{n+1},hf_{n}\big]^\top. \end{align*} $

将向前 Euler 方法(4.3)嵌入到上述 TSRK 中得到如下松弛 TSRK 方法 (RTSRK):

$\begin{matrix} {y}_{n+1}=y^{[n+1]}_1+\lambda_n(\tilde{y}^{[n]}_1-y^{[n+1]}_1). \end{matrix}$

参照文献[19], 离心率取值为 $e=0.6$, 分别用 RTSRK 和 TSRK 这两种方法在 $[0,100]$ 的时间区间内求解方程(4.7), 相应的数值结果详见表2. 从表2可以看出用 TSRK 方法进行长时间求解时, 大步长不能保证收敛阶为三阶, 而 RTSRK 方法能够适用更大的步长, 且在相同步长情形, 其精度要高于 TSRK, 且保持三阶收敛. 由于 TSRK 是三阶的, 向前 Euler 是一阶的, 所以 $\lambda_n$ 是关于 $h$ 的二阶无穷小量. 表2右半部分的数值结果验证了该理论结果的正确性. 同样, 我们也测试了方法的保结构性能, 详见图3. 从图3可知 RTSRK 相比于 TSRK 能更加精确地保持 Hamilton 量的不变性. 此外, 我们还绘制了$h=1/20$ 计算得到的数值解轨迹图, 详见图4. 从图4可知, 基于 TSRK 方法的数值解偏离了正确的轨道, 而基于 RTSRK 方法的数值解能够始终位于正确的椭圆轨道上.

表2   TSRK 方法和 RTSRK 方法对应的误差和收敛阶

新窗口打开| 下载CSV


图3

新窗口打开| 下载原图ZIP| 生成PPT

图3   $h=1/20$ 时 Kepler 方程的 Hamilton 量的相对误差演化图


图4

新窗口打开| 下载原图ZIP| 生成PPT

图4   解的轨迹图


例 4.3 考虑 sine-Gordon 方程:

$\begin{matrix} \begin{cases} &u_{tt}-u_{xx}+\sin u=0,~~(x,t)\in [-L,L]\times [T]\\ &u(x,0)=0,\\ &u_t(x,0)=4\kappa\mathrm{sech}(\kappa x),\\ &u(L,t)=u(-L,t), \end{cases} \end{matrix}$

其中 $\kappa=\frac{1}{\sqrt{1+c^2}}$, 相应的 Hamilton 量

$H(u(t))=\frac{1}{2}\int^L_{-L}|u_t|^2+|u_x|^2+2(1-\cos(u))\mathrm{d}x.$

根据文献[21]知, 该方程的精确解为:

$u(x,t)=4\arctan(c^{-1}\sin(c\kappa t)\mathrm{sech}(\kappa x)).$

在下面的计算过程中, 取参数 $c=0.5,~L=20$ 和 $T=100$.

令 $v=u_t$, sine-Gordon 方程(4.10)可改写为

$ \begin{align*} \begin{cases} &u_t=v, \\ &v_t=u_{xx}-\sin u. \end{cases} \end{align*} $

在空间上采取均匀网格剖分. 记步长 $h=2L/N$, 以及

$U(t)=\left[u_1,u_2,\cdots,u_N\right]^\top, V(t)=\left[v_1,v_2,\cdots,v_N\right]^\top.$

其中 $u_i\approx u(x_i,t),v_i\approx v(x_i,t)$ 和 $x_i=-L+(i-1)h$.空间二阶导数可以用二阶拟谱微分矩阵 $\mathcal{D}$ 逼近, 即

$\mathcal{D}=\mathcal{F}^H\Lambda \mathcal{F},$

其中 $\mathcal{F}$ 是离散 Fourier 矩阵, $\mathcal{F}^H$ 是 $\mathcal{F}$ 的共轭转置,

$\Lambda=-(\frac{2\pi}{2L})^2\text{diag}\left[0^2,1^2,\cdots,(\frac{N}{2})^2,(-\frac{N}{2}+1)^2,\cdots,(-2)^2,(-1)^2\right].$

关于谱方法更详细的介绍参考文献[22]. 相应的空间半离散系统可表述为

$\begin{matrix} \begin{cases} &U_t(t)=V(t), \\ &V_t(t)=\mathcal{D}U(t)-\sin(U(t)). \end{cases} \end{matrix}$

Hamilton 量的半离散形式为

$E(U(t),V(t))=\frac{h}{2}\left(V(t)^\top V(t)-U(t)^\top \mathcal{D}U(t)+2\textbf{e}^\top(\textbf{e}-\cos(U(t))) \right),$

其中 $\textbf{e}=\left[1,1,\cdots,1\right]^\top\in \mathbb{R}^N$.

考虑三阶显式 RK 方法:

$\begin{equation} \begin{cases}\label{RK3} Y_1=y^{[n]}_1+\frac12y^{[n]}_4,\\ Y_2=y^{[n]}_1-y^{[n]}_4+2\tau f(Y_1),\\ Y_3=\frac23\tau f(Y_1)+\frac16\tau f(Y_2)+y^{[n]}_1+\frac16y^{[n]}_4,\\ y^{[n+1]}_1=\frac23\tau f(Y_1)+\frac16\tau f(Y_2)+y^{[n]}_1+\frac16y^{[n]}_4,\\ y^{[n+1]}_i=y^{[n]}_{i-1}~~i=2,3,\\ y^{[n+1]}_4=\tau f(Y_3),\\ y^{[n+1]}_i=y^{[n]}_{i-1}.~~i=5,6, \end{cases} \end{equation}$

其中

$ \begin{align*} y^{[n]}=\big[y_{n},y_{n-1},y_{n-2},\tau f_{n},\tau f_{n-1},\tau f_{n-2}\big]^\top,~~ y^{[n+1]}=\big[y_{n+1},y_{n},y_{n-1},\tau f_{n+1},\tau f_{n},\tau y_{n-1}\big]^\top. \end{align*} $

考虑将三阶 Adams 外插法

$\begin{matrix} \hat{y}^{[n+1]}_1={y}^{[n]}_1+\frac{1}{12}(5{y}^{[n]}_4+8{y}^{[n]}_5-{y}^{[n]}_6)\end{matrix}$

嵌入到三阶 RK 方法(4.12)得到的松弛一般线性方法 (记为 RK-Adams 方法)

$\begin{matrix} y_{n+1}=y^{[n+1]}_1+\lambda_n(\hat{y}^{[n+1]}_1-y^{[n+1]}_1). \end{matrix}$

取空间步长 $h=2L/256$, 时间步长为 $\tau$, 分别应用该 RK 方法和 RK-Adams 方法去计算问题(4.11), 相应的计算结果详见表3和图5.

表3   RK-Adams 方法和 RK 方法的误差和收敛阶

新窗口打开| 下载CSV


图5

新窗口打开| 下载原图ZIP| 生成PPT

图5   $\tau=1/12$ 时 sine-Gordon 方程对应 Hamilton 量的相对误差演化图


从表3可知, RK-Adams 方法能够大幅提高计算精度. 从图5可知, 相较于单纯的 RK 方法, 松弛的一般线性方法 (RK-Adams 方法) 能在长时间模拟过程中更加精确地保持 Hamilton 量的不变性. 这再次验证了前面理论结果的正确性和本文所提出方法的有效性.

此外, 我们也将 RK-Adams 方法与二阶平均向量场方法 (AVF)[23] 以及 Gonzalez 离散梯度方法[8]进行了比较, 详见图6. 容易观察到, 与 AVF 方法和离散梯度方法相比, 基于本文思想构造的 RK-Adams 方法能够在相对非常短的计算时间内将误差降至更低的数量级. 这意味着, RK-Adams 方法在计算效率方面也显著优于 AVF 方法和离散梯度方法, 体现了在高维问题中的适用性.

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6   $T=10$ 时, 不同方法的计算效率对比图


5 总结

本文提出了一类基于松弛技术的线性隐式保结构算法, 并从理论层面对算法的误差和稳定性进行了分析. 结果表明, 该算法不仅能够精确守恒系统的不变量, 也在计算效率上具有显著优势, 且还保持了原方法的绝对稳定域. 最后, 基于三个代表性模型的数值结果也验证了理论分析结果的正确性以及算法的有效性. 需要指出的是, 本文的理论分析主要针对单不变量情形, 对多不变量系统的同步保持问题尚未深入探讨; 同时, 尽管数值实验验证了方法的稳定性, 其在刚性问题中的适用性仍有待进一步研究. 未来工作将致力于将该框架推广至多不变量系统, 并结合自适应步长策略以提升算法在非均匀动力学行为下的计算效率.

参考文献

Cooper G J.

Stability of Runge-Kutta methods for trajectory problems

IMA Journal of Numerical Analysis, 1987, 7(1): 1-13

DOI:10.1093/imanum/7.1.1      URL     [本文引用: 1]

Hairer E, Lubich C, Wanner G.

Geometric Numerical Integration:Structure-Preserving Algorithms for Ordinary Differential Equations

Berlin: Springer, 2006

[本文引用: 1]

Feng K, Qin M.

Symplectic Geometric Algorithms for Hamiltonian Systems

Berlin: Springer, 2010

[本文引用: 1]

Ruth R D.

A canonical integration technique

IEEE Transactions on Nuclear Science, 1983, 30: 2669-2671

DOI:10.1109/TNS.1983.4332919      URL     [本文引用: 1]

Marsden J E, West M.

Discrete mechanics and variational integrators

Acta Numerica, 2001, 10: 357-514

DOI:10.1017/S096249290100006X      URL     [本文引用: 1]

Brugnano L, Caccia G F, Iavernaro F.

Energy conservation issues in the numerical solution of the semilinear wave equation

Applied Mathematics and Computation, 2015, 270: 842-870

DOI:10.1016/j.amc.2015.08.078      URL     [本文引用: 1]

Simo J C, Tarnow N.

The discrete energy-momentum method

Zeitschrift Für Angewandte Mathematik Und Physik Zamp, 1992, 43: 757-792

[本文引用: 1]

Gonzalez O.

Time integration and discrete Hamiltonian systems

Journal of Nonlinear Science, 1996, 6: 449-467

DOI:10.1007/BF02440162      URL     [本文引用: 2]

McLachlan R I, Quispel G R W, Robidoux N.

Geometric integration using discrete gradients

Philos Trans A Math Phys Eng Sci, 1999, 357(1754): 1021-1045

DOI:10.1098/rsta.1999.0363      URL     [本文引用: 1]

Shampine L F.

Conservation laws and the numerical solution of ODEs

Computers & Mathematics with Applications, 1986, 12(5-6): 1287-1296

DOI:10.1016/0898-1221(86)90253-1      URL     [本文引用: 1]

Gear C W.

Maintaining solution invariants in the numerical solution of ODEs

SIAM Journal on Scientific and Statistical Computing, 1986, 7(3): 734-743

DOI:10.1137/0907050      URL     [本文引用: 1]

Del Buono N, Mastroserio C.

Explicit methods based on a class of four stage fourth order Runge-Kutta methods for preserving quadratic laws

Journal of Computational and Applied Mathematics, 2002, 140(1-2): 231-243

DOI:10.1016/S0377-0427(01)00398-3      URL     [本文引用: 1]

Ketcheson D I.

Relaxation Runge-Kutta methods: Conservation and stability for inner-product norms

SIAM Journal on Numerical Analysis, 2019, 57(6): 2850-2870

DOI:10.1137/19M1263662      [本文引用: 2]

We further develop a simple modification of Runge-Kutta methods that guarantees conservation or stability with respect to any inner-product norm. The modified methods can be explicit and retain the accuracy and stability properties of the unmodified Runge-Kutta method. We study the properties of the modified methods and show their effectiveness through numerical examples, including application to entropy-stability for first-order hyperbolic PDEs.

Li D, Li X, Zhang Z.

Implicit-explicit relaxation Runge-Kutta methods: Construction, analysis and applications to PDEs

Mathematics of Computation, 2023, 92: 117-146

DOI:10.1090/mcom/2023-92-339      [本文引用: 1]

谷伟, 李丁方, 李晓西, 张智民.

松弛隐显 Runge-Kutta 方法及其在高振荡 Hamilton 系统的应用

中国科学: 数学, 2025, 55(4): 829-848

[本文引用: 1]

Gu W, Li D F, Li X X, Zhang Z M.

Relaxation implicit-explicit Runge-Kutta method and its applications in highly oscillatory Hamiltonian systems (in Chinese)

Scientia Sinica Mathematica, 2025, 55(4): 829-848

DOI:10.1360/SSM-2023-0157      URL     [本文引用: 1]

Calvo M, Hernández-Abreu D, Montijano J I, et al.

On the preservation of invariants by explicit Runge-Kutta methods

SIAM Journal on Scientific Computing, 2006, 28(3): 868-885

DOI:10.1137/04061979X      URL     [本文引用: 3]

Butcher J C. Numerical Methods for Ordinary Differential Equations. New Jersey: John Wiley & Sons, 2016

[本文引用: 3]

Butcher J C, O'Sullivan A E.

Nordsieck methods with an off-step point

Numerical Algorithms, 2002, 31: 87-101

DOI:10.1023/A:1021104222126      [本文引用: 1]

Lu N, Cai W, Bo Y, et al.

Superconvergence of projection integrators for conservative system

Journal of Computational Physics, 2023, 490: Art 112281

[本文引用: 2]

Jackiewicz Z, Renaut R A, Zennaro M.

Explicit two-step Runge-Kutta methods

Applications of Mathematics, 1995, 40(6): 433-456

DOI:10.21136/AM      URL     [本文引用: 1]

Bratsos A G.

A numerical method for the one-dimensional sine-Gordon equation

Numerical Methods for Partial Differential Equations, 2008, 24(3): 833-844

DOI:10.1002/num.v24:3      URL     [本文引用: 1]

Shen J, Tang T, Wang L L.

Spectral Methods: Algorithms, Analysis and Applications

Berlin: Springer Science & Business Media, 2011

[本文引用: 1]

Celledoni E, Grimm V, McLachlan R I, et al.

Preserving energy resp. dissipation in numerical PDEs using the "Average Vector Field" method

Journal of Computational Physics, 2012, 231(20): 6770-6789

DOI:10.1016/j.jcp.2012.06.022      URL     [本文引用: 1]

/

〈 〉

版权所有 © 《数学物理学报》杂志社
地址：武汉市71010号信箱(430071)　电话: 027-87199206(中文版)　027-87199087(英文版)
E-mail: actams@wipm.ac.cn
本系统由北京玛格泰克科技发展有限公司设计开发