单延迟分段连续微分方程的数值稳定性

骆志纬

（广东工业大学应用数学学院，广东广州510520）

分段连续型微分方程（EPCAs）的基本形式为

这里 h（t）是常数区间，例如 h（t）=［t］，h（t）=［t-n］，h（t）=t-n［t］，其中 n 为正整数，［t］是不超过 t的最大整数。Cooke等［1］给出了该方程在生物模型和控制系统中的应用；Wiener［2］给出了关于该方程的一般理论与基础的结果，如该方程解的存在唯一性、稳定性、振动性和周期性。

2004年，刘明珠等［3］获得了如下分段连续型微分方程

的Runge-Kutta方法的数值解的稳定性结论，并应用（r，s）-Pade′逼近及Order-Star理论分析了Runge-Kutta方法的渐近稳定性。

2005年，宋明辉等［4］给出了如下分段连续型微分方程

的Runge-Kutta方法的数值稳定性。确定了Runge-Kutta方法的稳定区域，进一步得到了数值解稳定区域包含解析解稳定区域的充要条件。

2007年，刘明珠等［5］给出了如下分段连续型微分方程

的θ-方法的数值稳定性，同时也得到了该方程的方程数值解的稳定区域包含解析解的稳定区域的充要条件。

2008年，张如［6］应用θ-方法研究如下分段连续型微分方程

给出了数值解稳定区域包含解析解稳定区域的一个充分必要条件。接着，又应用线性θ-方法求解了方程x′（t）=au（t）+a1x（［t+p］）（p为正整数），给出了此类方程数值方法是渐近稳定的一个充分条件，同时也得出了数值解的稳定区域包含解析解的稳定区域的充分条件。

2016年，谢钰程等［7］研究了如下多维分段连续型微分方程

其中，L和M是n×n矩阵，利用无界区域最大模原理证明其差分格式为舒尔多项式，得到了θ-方法的稳定性结论。更多结果见文献［8-12］。

本文主要应用Runge-Kutta方法研究以下分段连续微分方程的数值稳定性

其中，a，a1，A0是实常数，［.］表示向下取整数函数。

定义 1［2］方程（1）的解在区间｛-1，0｝∪（0，∞）满足以下条件：

（1）u′（t）在［0，∞）上连续；

（2）导数u′（t）在每个点t∈［0，+∞）都存在，除了在点t∈［0，+∞）只有单侧导数；

（3）方程（1）在每个区间［n，n+1）∈［0，+∞）和整数端点成立。

定理 1［2］对任何给定的 A0和 A1，方程（1）在区间［0，+∞）存在唯一解

定理 2［3］对任意的实数 a，当

方程（4）全部根的模小于1。

定理3［2，6］方程（1）的解析解是渐近稳定的充分条件是

1 Runge-Kutta方法的稳定性

在这一节，考虑Runge-Kutta方法的稳定性。取步长为h=1/m，整数m≥1，网格点tn=nh（n=0，1，2，…）。假设b1+b2+…+bv=1且0≤c1≤c2≤…≤cv≤1，精确解u（t）在网格点tn（n=1，2，3，…）时的近似值为u1，u2，u3，…，那么Runge-Kutta方法对方程（1）的逼近可以得出以下形式的离散过程

其中，un是u（t）在tn（n=0，1，2，…）时的一个逼近值，yi（n）和zi（n）分别是u（tn+cih）和u（［tn+cih］）的逼近值。定义n=km+l（l=0，1，…，m-1）可以被定义为ukm，i=1，2，…，v，l=0，1，…，m-1。定义方程（7）归纳为

其中，l=0，1，…，m-1，e=（1，1，…，1）T。因此，当 a=0 时，有

当 a≠0 时，有

其中，l=0，1，…，m-2，x=ha，R（x）=1+xb（TI-xA）-1e是Runge-Kutta方法的稳定函数。设 φ（x）=b（TI-xA）-1e，因为 φ（0）=1 且 φ（x）是连续的，因此存在 δ＞0 使得对于所有x ≤δ，I-xA 是可逆的且 φ（x）＞0。

（1）（I-xA）是可逆的；

（2）对于任意给定的ui（1≤i≤m）由关系式（9）和（10）所确定的uk（k=1，2，…）满足当k→∞时，uk→0。

对于方程（1），使式（7）渐近稳定的所有点（a，a1）的集合称为数值解渐近稳定区域，记为S，把满足关系式（6）的解析解渐近稳定性的所有点（a，a1）组成的集合记为H，把集合H分为以下3个部分

当 a=0 时，由式（9）得到|ukm+l+1|≤|ukm|+|（l+1）ha1u（k+3）m|≤|ukm|+|a1||u（k+3）m|，故有 n→∞ 时，un→∞ 当且仅当 k→∞ 时 ukm→0。又 u（k+1）m=ukm+a1u（k+3）m，其特征方程与式（4）相同，故有 H⊂S。

故研究差分方程

该方程的特征方程为

要使ukm→0（k→∞）当且仅当方程（11）的所有特征根λ满足λ ＜1。

定理4 方程（1）的Runge-Kutta方法是渐近稳定的，当

令

记满足关系式（12）的所有（a，a1）的集合为 S1。

引理 1［14-17］ez的（r，s）-Pade′逼近由下式给出

其中

误差为

它是ez的唯一的阶为r+s的分子和分母的次数分别为r和s的有理逼近。

由文献［13-16］定义阶星

引理2［13-16］如果Runge-Kutta方法的阶为p，则当z→0时，D的形状类似于星型，由宽度都是πp+1的p+1个星形部分组成，这些星形部分被具有相同宽度的p+1个白色部分隔开。

引理 3［13-16］如果 R（z）是 ez的（r，s）-Pade′逼近，则：

（1）在右半平面有s个星形部分，每个都包含R（z）的一个极点；

（2）在左半平面有r个星形部分，每个都包含R（z）的一个零点；

（3）所有的部分都关于实轴对称。

推论1 假设R（z）是ez的（r，s）-Pade′逼近，如果s是奇数，那么在右半平面正实轴包含在星形部分中，由此得出对于所有的x＞0，R（x）＞ex；如果r是奇数，那么在左半平面负实轴包含在白色部分中，由此得出对于所有的 x＜0，R（x）＜ex。

定理5 假设方程（1）的Runge-Kutta方法的稳定函数由ex的（r，s）-pade逼近给出，那么当且仅当a=0 时，对任意的 r，s，有 H0⊂S；当 a＞0，s为奇数时，有 H1⊂S；当 a＜0，r为奇数时，有 H2⊂S。

（1）当a=0时，由定理3有

条件满足定理4，由推论1，故H0⊂S。

（2）当a＞0时，得Rm（x）≥ea＞1，有

故由定理3当

时，条件满足定理4，H1⊂S，从而由推论1即得s为奇数。

（3）当a＜0时，得1＞ea≥Rm（x）＞0，由上述（2）同理可得

时，条件满足定理4，H2⊂S，从而由推论即得r为奇数。

2 数值实验

给出如下数值例子，考虑方程

其中，a=2，a1=-3。

对于m=100，有R（x）=1+xbT（I-xA）-1e=1.020 201 340 026 756，且

对于m=50，有R（x）=1+xbT（I-xA）-1e=1.040 810 774 192 390，且

满足定理4，图1和图2画出了h=1/100和h=1/50时的Gauss-Legendre方法的数值解的图像，图像也说明了a与a1在满足定理4的条件下，数值解是渐近稳定的.

图1 h=1/100时方程（15）Gauss-Legendre方法的数值解

图2 h=1/50时方程（15）Gauss-Legendre方法的数值解

考虑如下方程

其中，a=-1，a1=3。

对于m=100，有R（x）=1+xbT（I-xA）-1e=0.990 049 833 749 168，且

对于m=50，有R（x）=1+xbT（I-xA）-1e=0.980 198 673 306 755，且

满足定理4，图3和图4画出了h=1/100和h=1/50时的LobattoⅢC方法的数值解的图像，图像也说明了a与a1在满足定理4的条件下，数值解是渐近稳定的。

图3 h=1/100时方程（16）LobattoⅢC方法的数值解

图4 h=1/50时方程（16）LobattoⅢC方法的数值解

考虑如下方程

图5 h=1/100时方程（17）RadauⅠA方法的数值解

图6 h=1/50时方程（17）RadauⅠA方法的数值解