随机延迟微分方程SK-ROCK方法的延迟相关稳定性分析

Delay Dependent Stability Analysis of SK-ROCK Methods for Stochastic Delay Differential Equations

摘要: 本文调整了显式稳定方法，即随机第二类正交Runge-Kutta-Chebyshev (SK-ROCK)方法，用于求解Itô随机延迟微分方程(SDDEs)。SK-ROCK方法实现了沿负实轴的扩展均方稳定区域，解决了现有显式方法在处理刚度和延迟相互作用方面的局限性。通过采用根定位技术，我们严格推导出了所提方法依赖于延迟的渐近均方稳定性条件。数值实验验证了SK-ROCK的收敛阶次和增强稳定性能，证明了它在实际应用中的有效性。

Abstract: This paper adapts the explicit stabilized method, termed the Stochastic Second kind Orthogonal Runge-Kutta-Chebyshev (SK-ROCK) method, for solving Itô stochastic delay differential equations (SDDEs). The SK-ROCK method achieves an extended mean-square stability region along the negative real axis, addressing the limitations of existing explicit methods in handling stiffness and delayed interactions. By employing root locus techniques, we rigorously derive the delay-dependent asymptotic mean-square stability conditions of the proposed method. Numerical experiments validate both the convergence order and enhanced stability performance of SK-ROCK, demonstrating its efficacy in practical applications.

1. 引言

随机延迟微分方程通过结合系统的随机扰动与历史依赖性，为复杂动态系统的建模提供了强有力的数学框架。作为随机微分方程与延迟微分方程的广义扩展，随机延迟微分方程在生物学、经济学、金融学及控制理论等领域具有广泛应用[1]-[3]。然而其解析解通常难以显式表达，使得数值方法的开发及其稳定性分析成为理论与应用研究的核心课题。

近年来，随机延迟微分方程数值解法的收敛性与稳定性研究取得了显著进展。Baker与Buckwar [4] [5]系统建立了显式一步法的强收敛性理论。Mao [6]在局部Lipschitz条件下证明了Euler-Maruyama方法的收敛性。Hu，Mohammed与Yan [7]提出了强一阶Milstein方法。在稳定性分析中，Huang等人[8]揭示了随机Theta方法的延迟相关渐进均方稳定性特性，而Hu与Huang [9]通过随机指数欧拉方法，推导了延迟相关稳定性的充要条件，深化了对数值方法鲁棒性的理解。

在刚性随机延迟微分方程中，数值求解常面临显式方法稳定性受限与隐式方法计算成本高昂的矛盾。Euler-Maruyama方法等经典显式方法需采用极小步长以保证稳定性，而半隐式或隐式方法虽具备优越稳定性，却因需迭代求解非线性方程组而难以适用于大规模系统。为此，研究者致力于发展兼顾效率与鲁棒性的显式算法。随机正交Runge-Kutta-Chebyshev (S-ROCK)方法[10]-[14]通过显式结构设计在随机常微分方程中实现了稳定性突破；Guo等[15]将其扩展至随机延迟微分方程并提出显式Runge-Kutta-Maruyama方法，Komori等人[16]进一步优化了固定延迟情形下的阻尼因子。然而，现有方法仍存在稳定域受限、延迟适应性不足等问题，尤其是经典S-ROCK方法中随阶段数递增的阻尼因子η虽抑制了噪声扰动，却显著缩小了稳定域长度。因此，设计具有扩展稳定域的高效显式方法成为当前研究的关键方向。Abdulle等人[17]提出了一种具有最优扩展均方稳定域的弱一阶显式稳定方法。

本文将一类显式随机第二类正交Runge-Kutta-Chebyshev (SK-ROCK)方法扩展到SDDEs。其目的是通过优化阻尼机制和延迟项处理策略，大幅扩展均方稳定域。文章第2节回顾切比雪夫方法的基本理论；第3节推导SK-ROCK方法的构造流程；第4节系统分析其延迟相关稳定性准则；第5节通过数值实验验证理论结果；第六节给出本文的结论。

2. 切比雪夫方法

考虑常微分方程初值问题：

$\text{d}y\left(t\right)=f\left(y\left(t\right)\right)\text{d}t,t\in \left[0,T\right],y\left(0\right)=y_0.$ (1)

van der Houwen和Sommeijer [18]提出了一类基于s级显式Runge-Kutta框架的稳定算法。该方法通过引入第一类切比雪夫多项式构造稳定性多项式，因此被称为切比雪夫方法。其核心定义如下：

$K_0=y_n,K_1=y_n+h\frac{{\omega }_1}{{\omega }_0}f\left(K_0\right),$

$K_i=2\frac{T_{i-1}\left({\omega }_0\right)}{T_i\left({\omega }_0\right)}\left(h{\omega }_{1}f\left(K_{i-1}\right)+{\omega }_0{K}_{i-1}\right)-\frac{T_{i-2}\left({\omega }_0\right)}{T_i\left({\omega }_0\right)}{K}_{i-2},i=2,3,\cdots ,s,$ (2)

$y_{n+1}=K_s$

其中参数${\omega }_0$ 与${\omega }_1$ 满足

${\omega }_0=1+\frac{\eta }{s^2},{\omega }_1=\frac{T_s\left({\omega }_0\right)}{{T^{\prime }}_s\left({\omega }_0\right)}$ (3)

$T_k\left(x\right)$ 是阶数为k的切比雪夫多项式，具有如下递推关系：

$T_0\left(x\right)=1,T_1\left(x\right)=x,T_k\left(x\right)=2x{T}_{k-1}\left(x\right)-T_{k-2}\left(x\right),k=2,3,\cdots ,s$ (4)

其中$y_n$ 表示$y\left(t\right)$ 在等距步长点$t_n=nh\left(0\le n\le N\right)$ 的数值近似值，并且对于给定的正整数$N$ ，时间步长等于$h=T/N$ 。

将方法(2)应用于测试方程

$y^{\prime }\left(t\right)=\lambda y\left(t\right),y\left(0\right)=y_0$ (5)

其中$\Re \left(\lambda \right)\le 0$ ，$y_0\ne 0$ ，可以得到

$y_{n+1}=R_s\left(\lambda h\right)y_n=\frac{T_s\left({\omega }_0+{\omega }_1\lambda h\right)}{T_s\left({\omega }_0\right)}{y}_n$ (6)

$R_s\left(z\right)$ 称为该数值方法的稳定函数。该方法的稳定域为$S=\left\{z:\left|R_s\left(z\right)\right|\le 1,z\in ℂ\right\}$ 。当阻尼值$\eta =0$ 时，稳定函数可简化为$R_s\left(z\right)=T_s\left(1+z/s^2\right)$ 。根据切比雪夫多项式的极值性质，当$z\in \left[-2s^2,0\right]$ 时，稳定函数满足$\left|R_s\left(z\right)\right|\le 1$ 。因此，随着$R_s\left(z\right)$ 的级数s的增加，该方法的稳定域沿着负实轴作二次扩展并满足

$\left[-l_s^{\eta }\cdot s^2,0\right]\subset S$

其中稳定域长度系数$l_s^{\eta }\cong 2-4/3\eta$ 。特别地，当$\eta =0$ 时，$l_s^{\eta }$ 达到最大值2，此时稳定域覆盖负实轴区间$\left[-2s^2,0\right]$ 。

由图1和图2可知，切比雪夫方法的稳定域由条件$\left|R_s\left(z\right)\right|\le 1$ 所定义，其在负实轴上的覆盖范围随阻尼参数$\eta$ 变化。当$\eta =0$ 时，稳定域沿负实轴达到最大区间$\left[-2s^2,0\right]$ ，但此时稳定域在虚轴附近退化为离散点，导致数值方法对虚轴方向的特征值扰动极其敏感。为此，引入适量阻尼(如$\eta =0.05$ )可增强鲁棒性，尽管稳定域总长度略微缩减，但其覆盖范围从纯负实轴向邻近复平面区域扩展，从而有效抑制因高频扰动引发的数值失稳现象。

Figure 1. Complex stability domains of Chebyshev method for $s=6$ with $\eta =0,0.05$

图1. 当$s=6$ ，$\eta =0,0.05$ 时切比雪夫方法的复稳定域

Figure 2. Stability function of Chebyshev method for $s=6$ with $\eta =0,0.05$

图2. 当$s=6$ ，$\eta =0,0.05$ 时切比雪夫方法的稳定函数

Abdulle和Li [11]进一步将切比雪夫方法推广至Itô型随机微分方程

$\text{d}y\left(t\right)=f\left(y\left(t\right)\right)\text{d}t+g\left(y\left(t\right)\right)\text{d}W\left(t\right),y\left(0\right)=y_0,$ (7)

并提出显式随机正交Runge-Kutta-Chebyshev (S-ROCK)方法：

$K_0=y_n,K_1=y_n+h\frac{{\omega }_1}{{\omega }_0}f\left(K_0\right),$

$K_i=2\frac{T_{i-1}\left({\omega }_0\right)}{T_i\left({\omega }_0\right)}\left(h{\omega }_{1}f\left(K_{i-1}\right)+{\omega }_0{K}_{i-1}\right)-\frac{T_{i-2}\left({\omega }_0\right)}{T_i\left({\omega }_0\right)}{K}_{i-2},i=2,3,\cdots ,s,$ (8)

$y_{n+1}=K_s+g\left(K_s\right)\Delta W$

在确定性方法(2)中，$\eta$ 被选择得很小且固定不变，通常取$\eta =0.05$ ，而在S-ROCK方法(8)中，阻尼通常取较大值并且是s的递增函数。为了获得最佳稳定域，Abdulle等人[17]提出了SK-ROCK方法：

$K_0=y_n,$

$K_1=y_n+h\frac{{\omega }_1}{{\omega }_0}f\left(K_0+\frac{s{\omega }_1}{2}g\left(K_0\right)\Delta W\right)+\frac{s{\omega }_1}{{\omega }_0}g\left(K_0\right)\Delta W,$

$K_i=2\frac{T_{i-1}\left({\omega }_0\right)}{T_i\left({\omega }_0\right)}\left(h{\omega }_{1}f\left(K_{i-1}\right)+{\omega }_0{K}_{i-1}\right)-\frac{T_{i-2}\left({\omega }_0\right)}{T_i\left({\omega }_0\right)}{K}_{i-2},i=2,3,\cdots ,s,$ (9)

$y_{n+1}=K_s.$

3. SDDE的SK-ROCK方法

设$\left(\Omega ,\mathfrak{F},{\left\{{\mathfrak{F}}_t\right\}}_{t\ge 0},\mathbb{P}\right)$ 是一个完备的概率空间，其中滤子${\left\{{\mathfrak{F}}_t\right\}}_{t\ge 0}$ 右连续且满足通常条件(即${\mathfrak{F}}_0$ 包含所有$\mathbb{P}$ -零集)。考虑以下形式的Itô型随机延迟微分方程：

$\{\begin{array}{l}\text{d}y\left(t\right)=f\left(y\left(t\right)\right)\text{d}t+g\left(y\left(t\right),y\left(t-\tau \right)\right)\text{d}W\left(t\right),t\in \left[0,T\right],\\ y\left(t\right)=\Psi \left(t\right),t\in \left[-\tau ,0\right]，\end{array}$ (10)

其中$\tau >0$ 为常数延迟，$f:{ℝ}^d\times {ℝ}^d\to {ℝ}^d$ 是漂移项，$g:{ℝ}^d\times {ℝ}^d\to {ℝ}^d$ 是扩散项，$W\left(t\right)$ 为定义在上述概率

空间的一维标准维纳过程。初始函数$\Psi \left(t\right)$ 在$\left[-\tau ,0\right]$ 上连续，并满足 $\text{E}\left[{\text{sup}}_{t\in \left[-\tau ,0\right]}{\left|\Psi \left(t\right)\right|}^2\right]<+\infty$ 。为了保证

方程(10)的解的存在性与唯一性，假定漂移项和扩散项足够光滑并且满足以下条件[5]：

1) (全局Lipschitz条件)存在常数$L_{1,2}>0$ ，使得$\forall t\in \left[0,T\right]$ 以及$\forall x,y,\overline{x},\overline{y}\in {ℝ}^d$ ，有

$\left|f\left(x\right)-f\left(\overline{x}\right)\right|\le L_1\left(\left|x-\overline{x}\right|\right),\left|g\left(x,y\right)-g\left(\overline{x},\overline{y}\right)\right|\le L_2\left(\left|x-\overline{x}\right|+\left|y-\overline{y}\right|\right)$ (11)

2) (线性增长条件)存在常数$K_{1,2}>0$ ，使得$\forall t\in \left[0,T\right]$ 以及$\forall x,y\in {ℝ}^d$ ，有

${\left|f\left(x\right)\right|}^2\le K_1\left(1+{\left|x\right|}^2\right),{\left|g\left(x,y\right)\right|}^2\le K_2\left(1+{\left|x\right|}^2+{\left|y\right|}^2\right)$ (12)

本文假设延迟$\tau =Mh$ ，其中时间步长$h=T/N$ (M、N均为正整数)。(10)的显式EM方法如下[19]：

$y_{n+1}=y_n+f\left(y_n\right)h+g\left(y_n,y_{n-M}\right)\Delta W_n.$ (13)

其中$\Delta W_n=W\left(t_{n+1}\right)-W\left(t_n\right)$ 。初始阶段($n=0,1,\cdots ,M$ )的延迟项$y_{n-M}$ 由初始函数$\Psi \left(t_n-\tau \right)$ 确定，后续阶段($n\ge M+1$ )则由数值解递推生成。将显式随机第二类正交Runge-Kutta-Chebyshev (SK-ROCK)方法扩展到(10)。那么SK-ROCK方法可以写成

$K_0=y_n,K_0^{\text{*}}=K_0+\frac{s{\omega }_1}{2}g\left(K_0,y_{n-M}\right)\Delta W_n,$

$K_1=y_n+h\frac{{\omega }_1}{{\omega }_0}f\left(K_0^{\text{*}}\right)+\frac{s{\omega }_1}{{\omega }_0}g\left(K_0,y_{n-M}\right)\Delta W_n,$

$K_i=2\frac{T_{i-1}\left({\omega }_0\right)}{T_i\left({\omega }_0\right)}\left(h{\omega }_{1}f\left(K_{i-1}\right)+{\omega }_0{K}_{i-1}\right)-\frac{T_{i-2}\left({\omega }_0\right)}{T_i\left({\omega }_0\right)}{K}_{i-2},i=2,3,\cdots ,s,$ (14)

$y_{n+1}=K_s,$

设$y\left(t\right)$ 为方程(10)的精确解，$y_n$ 为其在$t_n$ 处的数值近似。方法(14)在$t_n$ 时逼近精确解 $y\left(t\right)$ 的全局误差定义为${\epsilon }_n=y\left(t_n\right)-y_n$ ，如果存在正常数$c$ 和${\delta }_0$ ，使得

$\underset{0\le n\le N}{\max }{\left(E{\left|{\epsilon }_n\right|}^2\right)}^{1/2}\le c{h}^q,\forall h\in \left(0,{\delta }_0\right)$ (15)

则称该方法的强均方收敛阶为q。

定理3.1. 假设方程(10)满足条件(11)与(12)。则数值方法(14)的强均方收敛阶为1/2。

Komori、Eremin和Burrage [16]已证明S-ROCK方法具有强1/2阶收敛性，类似结论可推广至SK-ROCK方法。为简洁起见，本文省略了定理证明的细节，收敛结果将在第5节中通过数值实验进行验证。

4. 延迟相关稳定性

本节分析SK-ROCK方法的延迟相关渐进均方稳定性。考虑以下线性随机延迟方程

$\{\begin{array}{l}\text{d}y\left(t\right)=\lambda y\left(t\right)\text{d}t+\left(\mu y\left(t\right)+\sigma y\left(t-\tau \right)\right)\text{d}W\left(t\right),t\in \left[0,T\right],\\ y\left(t\right)=\Psi \left(t\right),t\in \left[-\tau ,0\right],\end{array}$ (16)

Appleby，Mao和Riedle [20]分析了方程(13)的渐近均方稳定性，得到以下命题。

命题4.1. 设$\lambda <0$ ，$\mu ,\sigma \in ℝ$ ，$\tau >0$ ，则方程(13)的渐近均方稳定区域$S_{*}$ 为

$S_{*}=\left\{\left(\lambda ,\mu ,\sigma \right):{\mu }^2+{\sigma }^2+2\mu \sigma {\text{e}}^{\lambda \tau }<-2\lambda \right\}$ (17)

将SK-ROCK方法(14)应用于上述方程(16)，得到以下递推形式

$y_{n+1}=A\left(p\right)y_n+B\left(p\right)\left(\mu y_n+\sigma y_{n-M}\right)\Delta W_n.$ (18)

其中系数函数为

$A\left(p\right)=\frac{T_s\left({\omega }_0+{\omega }_{1}p\right)}{T_s\left({\omega }_0\right)},B\left(p\right)=\frac{U_{s-1}\left({\omega }_0+{\omega }_{1}p\right)}{U_{s-1}\left({\omega }_0\right)}\left(1+\frac{{\omega }_1}{2}p\right).$

此处$T_s\left(x\right)$ 为第一类切比雪夫多项式，$U_s\left(x\right)$ 为第二类切比雪夫多项式，满足递推关系：

$U_0\left(x\right)=1,U_1\left(x\right)=2x,U_k\left(x\right)=2x{U}_{k-1}\left(x\right)-U_{k-2}\left(x\right),k=2,3,\cdots ,s$ (19)

对式(18)两边平方并取数学期望，可得

$E\left[y_{n+1}^2\right]=A{\left(p\right)}^{2}E\left[y_n^2\right]+B{\left(p\right)}^2\left(h{\mu }^{2}E\left[y_n^2\right]+h{\sigma }^{2}E\left[y_{n-M}^2\right]+2h\mu \sigma E\left[y_n{y}_{n-M}\right]\right),$ (20)

由于$\Delta W_n$ 与$y_n$ 和$y_{n-M}$ 独立，交叉项期望可递推化简为

$E\left[y_n{y}_{n-M}\right]=A{\left(p\right)}^{M}E\left[y_{n-M}^2\right],$

令$p=\lambda h$ ，$q=\mu \sqrt{h}$ ，$v=\sigma \sqrt{h}$ ，由(20)可得

$E\left[y_{n+1}^2\right]=A{\left(p\right)}^{2}E\left[y_n^2\right]+B{\left(p\right)}^2\left(h{\mu }^{2}E\left[y_n^2\right]+h{\sigma }^{2}E\left[y_{n-M}^2\right]+2h\mu \sigma A{\left(p\right)}^{M}E\left[y_{n-M}^2\right]\right).$ (21)

其相应的特征方程为

$\xi =A{\left(p\right)}^2+q^{2}B{\left(p\right)}^2+B{\left(p\right)}^2\left(v^2+2qvA{\left(p\right)}^M\right){\xi }^{-M}.$ (22)

当p、q、v和M给定时，如果(22)的所有根都满足$\left|\xi \right|<1$ ，则方法(14)具有渐近均方稳定性。因此，对于任意固定的M，方法(14)的均方稳定区域为

$S_M=\left\{\left(p,q,v\right):\left(22\right)的所有根\xi 满足\left|\xi \right|<1\right\}$

当$q=v=0$ 时，特征方程退化为$\xi =A{\left(p\right)}^2$ ，此时稳定性条件简化为$\left|A\left(p\right)\right|<1$ 。

基于多项式根连续依赖于系数的特性，采用根轨迹技术分析方程(22)的根分布确定参数平面$\left(q,v\right)$ 在固定$p$ 值下的稳定边界。首先分析模为1的特征根位置。设$\xi ={\text{e}}^{i\phi }$ ，$\phi \in ℝ$ ，由于特征方程(22)的系数都是实数，只需要考虑$\phi \in \left[0,\pi \right]$ 的情形。通过遍历$\phi$ 值域，可解析参数平面内的稳定边界曲线。

当$\phi =0$ 时，即$\xi =1$ ，由方程(22)可得曲线

$C_0=\left\{\left(q,v\right):q^2+v^2+2qvA{\left(p\right)}^M+\frac{A{\left(p\right)}^2-1}{B{\left(p\right)}^2}=0\right\}$ (23)

对于给定的$p$ 值，它表示$\left(q,v\right)$ 平面上的一条椭圆曲线。

当$\phi =\pi$ 时，即$\xi =-1$ ，需区分延迟步数M的奇偶性。当M为偶数时，方程(22)右侧可展开为

$A{\left(p\right)}^2+B{\left(p\right)}^2{\left(q+vA{\left(p\right)}^M\right)}^2+v^2\left(1-A{\left(p\right)}^{2M}\right),$

由于$A{\left(p\right)}^{2M}<1$ (由$\left|A\left(p\right)\right|<1$ 保证)，右侧表达式恒为正，故$\xi =-1$ 非方程根。当M为奇数时，临界曲线为双曲线族

$C_{\pi }=\left\{\left(q,v\right):q^2-v^2-2qvA{\left(p\right)}^M+\frac{A{\left(p\right)}^2+1}{B{\left(p\right)}^2}=0\right\}$ (24)

其分支$C_{\pi }^1$ (上半平面)与$C_{\pi }^2$ (下半平面)关于原点对称。

当$\phi \in \left(0,\pi \right)$ 时，将$\xi ={\text{e}}^{i\phi }$ 代入特征方程(22)，分离实部和虚部得到

$\cos \phi =A{\left(p\right)}^2+q^{2}B{\left(p\right)}^2+B{\left(p\right)}^2\left(v^2+2qvA{\left(p\right)}^M\right)\cos M\phi ,$ (25)

$\sin \phi =B{\left(p\right)}^2\left(v^2+2qvA{\left(p\right)}^M\right)\sin M\phi .$ (26)

由(26)可知，当$\xi ={\text{e}}^{ik\pi /M},k=1,\cdots ,M-1$ 时，$\sin M\phi =0$ 但$\sin \phi \ne 0$ ，此类根不存在。将区间$\phi \in \left(0,\pi \right)$

划分为M个子区间$\left(\frac{k\pi }{M},\frac{\left(k+1\right)\pi }{M}\right),k=0,1,\cdots ,M-1$ ，由(25)~(26)得

$q^2=\frac{\sin \left(M+1\right)\phi }{\sin M\phi }\cdot \frac{1}{B{\left(p\right)}^2}-\frac{A{\left(p\right)}^2}{B{\left(p\right)}^2},$ (27)

$v^2+2qvA{\left(p\right)}^M+\frac{\sin \phi }{\sin M\phi }\cdot \frac{1}{B{\left(p\right)}^2}=0.$ (28)

将这些子区间分为两类进行分析。首先考虑k为偶数的情况，这时$\sin M\phi >0$ ，当$\sin \left(M+1\right)\phi \le 0$ 时，式(27)无实数解q。当$\sin \left(M+1\right)\phi >0$ 时，方程(28)关于变量v的判别式$\text{Δ}$ 为

$\text{Δ}=4\frac{\sin \phi }{\sin M\phi }\cdot \frac{1}{B{\left(p\right)}^2}\left[A{\left(p\right)}^{2M}\frac{\sin \left(M+1\right)\phi }{\sin M\phi }-A{\left(p\right)}^{2M+2}\frac{\sin M\phi }{\sin \phi }-1\right],$ (29)

可以推断出判别式$\text{Δ}<0$ ，故无实数解v。因此，当k为偶数时，对于任意$\phi \in \left(\frac{k\pi }{M},\frac{\left(k+1\right)\pi }{M}\right)$ ，(27)~(28)关于q、v都没有实数解。

其次，考虑k为奇数的情况，此时$\sin M\phi <0$ 。由于函数$\frac{\sin \left(M+1\right)\phi }{\sin M\phi }$ 在区间$\left(\frac{k\pi }{M},\frac{\left(k+1\right)\pi }{M}\right)$ 内严格单调递减，则对于任何固定的M > 1和k，都存在一个唯一的临界角${\phi }_k \in \left(0,\pi \right)$ ，使得

$\frac{\sin \left(\left(M+1\right)\cdot \frac{k\pi +{\phi }_k}{M}\right)}{\sin \left(M\cdot \frac{k\pi +{\phi }_k}{M}\right)}=A{\left(p\right)}^2$

当$\phi \in \left(\frac{k\pi +{\phi }_k}{M},\frac{\left(k+1\right)\pi }{M}\right)$ 时，由${\phi }_k$ 的定义可知$\frac{\sin \left(M+1\right)\phi }{\sin M\phi }<A{\left(p\right)}^2$ ，因此(27)关于q无实数解。当$\phi \in \left(\frac{k\pi }{M},\frac{k\pi +{\phi }_k}{M}\right]$ 时，方程(27)存在一对实数解q，且二次方程(28)的判别式$\text{Δ}\ge 0$ ，故对应实解

$v=-qA{\left(p\right)}^M\pm \sqrt{\text{Δ}/4}$ 。由此生成上下半平面对称曲线

$C_k=\{\left(q,v\right):q^2=\frac{\sin \left(M+1\right)\phi }{\sin M\phi }\cdot \frac{1}{B{\left(p\right)}^2}-\frac{A{\left(p\right)}^2}{B{\left(p\right)}^2},v=-qA{\left(p\right)}^M+\sqrt{\frac{\text{Δ}}{4}},\phi \in (\frac{k\pi }{M},\frac{k\pi +{\phi }_k}{M}]\}$

$C_{k+1}=\{\left(q,v\right):q^2=\frac{\sin \left(M+1\right)\phi }{\sin M\phi }\cdot \frac{1}{B{\left(p\right)}^2}-\frac{A{\left(p\right)}^2}{B{\left(p\right)}^2},v=-qA{\left(p\right)}^M-\sqrt{\frac{\text{Δ}}{4}},\phi \in (\frac{k\pi }{M},\frac{k\pi +{\phi }_k}{M}]\}$

显然，$C_k$ 位于上半平面，$C_{k+1}$ 位于下半平面，并且相对于原点是相互对称的。[8]中给出的引理分析了上述曲线的性质。

引理4.2. (曲线相交性)曲线$C_k\left(k\ge 1\right)$ 互不相交，每条$C_k$ 与直线

$q=\sqrt{\frac{1}{B{\left(p\right)}^2}-\frac{A{\left(p\right)}^2}{B{\left(p\right)}^2}}$

相交一次，交点的$v$ 坐标$v_k$ 满足 $\left|v_{k+2}\right|>\left|v_k\right|$ ，表明交点模值随k递增。

引理4.3. (区域分布)曲线$C_0$ 与其余曲线不相交，且与q轴相交两次，位于$C_1$ 和$C_2$ 之间。当M为奇数时，曲线$C_{\pi }^1$ 位于$C_{M-2}$ 的上方，而$C_{\pi }^2$ 位于$C_{M-1}$ 的下方。

由引理4.2和引理4.3可知，曲线族$C_0$ 、$C_k$ 、$C_{\pi }^1$ 和$C_{\pi }^2$ 将参数平面$\left(q,v\right)$ 划分为$2\lfloor \left(M+3\right)/2\rfloor$ 个区域($\lfloor x\rfloor$ 为向下取整函数)。根据特征方程根的连续性，每个区域中单位盘外的根数是恒定的。显然，当$\left(q,v\right)=\left(0,0\right)$ 时，方程(22)的所有M + 1个根均位于单位圆内，故$C_0$ 界内区域为稳定域。此外，如果向上穿过曲线$C_k$ ，一对共轭复根进入右半平面，单位圆外根数增加2。同样地，如果向下穿过曲线$C_{k+1}$ ，单位圆外根数也增加2。

在$C_0$ 区域之外，$C_1$ 和$C_2$ 所限定的区域内，方程(22)仅有一个根的模大于1。此外，当m为奇数时，$C_{\pi }^1$ 和$C_{\pi }^2$ 存在，当点$\left(q,v\right)$ 位于曲线$C_{\pi }^1$ 的上方或曲线$C_{\pi }^2$ 的下方区域时，(22)的所有M + 1根都在单位圆之外。根据上述分析，SK-ROCK法的渐进均方稳定性可归纳为以下命题。

命题4.4. 如果$\lambda$ ，$\mu$ ，$\sigma$ 和$h$ 满足

$\left|A\left(p\right)\right|<1$ , $q^2+v^2+2qvA{\left(p\right)}^M+\frac{A{\left(p\right)}^2-1}{B{\left(p\right)}^2}<0$

则(16)的SK-ROCK方法是渐进均方稳定的。

若只考虑延迟随机项[16]，即$\mu =0$ ，方程(16)的精确解是均方稳定的当系数满足${\sigma }^2<-2\lambda$ 。此时，SK-ROCK方法的均方稳定域为

$S_M=\left\{\left(p,v\right):A{\left(p\right)}^2+B{\left(p\right)}^2{v}^2<1\right\}$

Figure 3. MS stability domains for SK-ROCK with $s=6$ when $\eta =0.05$ and 3.61

图3. 当$\eta =0.05$ 和3.61时，$s=6$ 的SK-ROCK方法的均方稳定域

Figure 4. MS stability domains for S-ROCK with $s=6$ when $\eta =0.05$ and 3.61

图4. 当$\eta =0.05$ 和3.61时，$s=6$ 的S-ROCK方法的均方稳定域

图3和图4分别给出了$s=6$ 时不同阻尼值下标准SK-ROCK方法和S-ROCK方法的稳定域，灰色部分表示数值方法的稳定域，网格围成的区域表示$\mu =0$ 时(16)的均方稳定域。与确定性的切比雪夫方法类似，假定存在一个值$-l_s^{\eta }\cdot s^2$ ，使得对于任意的$p>-l_s^{\eta }\cdot s^2$ ，数值方法的稳定域包含方程(16)的稳定域。如图4所示，标准S-ROCK方法存在一个关键的稳定性限制：较小的阻尼值(如$\eta =0.05$ )会产生违反准则的不稳定区域。尽管[16]通过增加阻尼参数解决了这一问题，例如在$s=6$ 时$\eta =3.61$ ，但如此高的阻尼大大减少了负实轴稳定域的长度。相比之下，图3说明了SK-ROCK方法的稳健性。在相同级数($s=6$ )和低阻尼值($\eta =0.05$ )条件下，SK-ROCK不仅能保持完全符合稳定性要求，还能在负实轴方向上大幅扩展稳定域。这些优势使它特别适用于求解具有快速衰减成分的刚性随机微分方程，而S-ROCK等传统方法很难同时保持稳定性和稳定域覆盖范围。

5. 数值实验

本节将给出一些数值实验来验证我们的理论结果。考虑以下线性随机延迟微分方程[21]

$\{\begin{array}{l}\text{d}y\left(t\right)=\left(ay\left(t\right)+by\left(t-1\right)\right)\text{d}t+\left(cy\left(t\right)+dy\left(t-1\right)\right)\text{d}W\left(t\right),t\ge 0,\\ y\left(t\right)=t+1,t\in \left[-1,0\right],\end{array}$ (30)

当$t\in \left[0,1\right]$ 时，(30)的解为

$y\left(t\right)={\varphi }_{t,0}\left(1+{\displaystyle {\int }_0^t{\varphi }_{s,0}^{-1}}\left(b-cd\right)s\text{d}s+{\displaystyle {\int }_0^{t}ds{\varphi }_{s,0}^{-1}\text{d}W\left(s\right)}\right),$ (31)

${\varphi }_{t,0}=\text{exp}\left({\displaystyle {\int }_0^t\left(a-\frac{1}{2}{c}^2\right)\text{d}s}+{\displaystyle {\int }_0^{t}c\text{d}W\left(s\right)}\right)$

当$t\in \left[1,2\right]$ 时，将(31)作为新的初始函数，从而得到显式解

$y\left(t\right)={\varphi }_{t,1}\left(y\left(1\right)+{\displaystyle {\int }_1^t{\varphi }_{s,0}^{-1}\left(b-cd\right)y\left(s-1\right)\text{d}s}+{\displaystyle {\int }_1^{t}d{\varphi }_{s,0}^{-1}y\left(s-1\right)\text{d}W\left(s\right)}\right),$

以同样的方法，可以一步一步地得到后续区间的显式解。

为验证SK-ROCK方法的理论收敛阶，取$a=-8$ ，$b=0$ ，$c=d=0.5$ ，在终态时间T = 2下，分别对EM方法及二阶SKROCK方法进行收敛性分析。均方根误差(RMSE)通过分块统计方法计算：将总样本划分为20个独立区块，每个区块包含100条样本路径，最终取各区块误差的均值以提升统计可靠性。图5展示了步长与RMSE的对数收敛关系。图中实线与虚线分别对应SK-ROCK方法与EM方法的实验结果，点线为斜率1/2的理论参考线。实验结果表明，两种方法的RMSE随步长减小的速率均与参考线吻合，验证了其收敛阶为理论预期的1/2阶。

与EM方法相比，SK-ROCK方法的优势在于稳定性。为了说明这一点，取$a=-10$ ，$b=0$ ，$c=d=1$ ，根据命题4.1的理论分析，该方程精确解具有渐近均方稳定性。通过10,000次蒙特卡洛模拟生成数值解样本路径，图6展示了不同步长下两种方法的均方范数演化规律。实验表明，当步长h = 1/4时，EM方法因稳定性限制导致均方范数发散，而SK-ROCK方法仍保持稳定衰减。

进一步考虑非线性随机延迟方程[8]：

$\text{d}y\left(t\right)=ay\left(t\right)\left(1-y\left(t-1\right)\right)\text{d}t+\left(cy\left(t\right)+dy\left(t-1\right)\right)\text{d}W\left(t\right),t\ge 0,$ (32)

Figure 5. RMSE of y(2). Solid line: SK-ROCK; dashed line: EM; dotted line: reference line with slope 1/2

图5. y(2)的均方根误差。实线：SK-ROCK；虚线：EM；点线：斜率为1/2的参考线

Figure 6. Numerical simulations with different step sizes. (left: EM method; right: SK-ROCK method)

图6. 不同步长下的数值模拟。(左：EM法；右：SK-ROCK法)

Figure 7. Numerical simulations with different step sizes. (Top Left: EM method; Top Right: S-ROCK method; Bottom: SK-ROCK method)

图7. 不同步长下的数值模拟。(左上：EM法；右上：S-ROCK法；底部：SK-ROCK法)

设$a=-15$ ，$c=d=1$ ，初始条件$y\left(t\right)=0.1,t\in \left[-1,0\right]$ 。我们在不同步长条件下，对10,000条数值模拟的样本路径进行均方平均分析，实验结果如图7所示：其中左侧图像表示采用EM方法的模拟结果；中间图像展示的是标准S-ROCK方法的模拟结果；而右侧图像则对应于SK-ROCK方法的模拟效果。实验结果显示，在所有测试的步长条件下，EM方法表现出不稳定性。同样地，尽管S-ROCK方法在较小步长(如1/4)下能保持稳定，但在更大步长下却出现了不稳定现象。相比之下，SK-ROCK方法无论在何种步长下均显示出了良好的稳定性。

6. 结论

根据本文的分析，SK-ROCK与S-ROCK相比，适用于扩散项有延迟的SDDE更稳定，所需的阻尼系数也更低，能提供良好的均方稳定性。这使它们成为SDEs和SDDEs的更好选择。不过，在这里我们有意忽略了漂移项中的延迟。在这种情况下，Runge-Kutta-Chebyshev方法的稳定性分析变得更加困难(见[22] [23])，这将是未来研究的主题。

基金项目

国家自然科学基金(12172186；11772166)。

NOTES

^*第一作者。

^#通讯作者。