中立型时滞微分方程Rosenbrock 方法的弱时滞相关稳定性

张明坤 王艳沛 赵欢欢

(1.上海大学理学院, 上海 200444;2.郑州航空工业管理学院数学学院, 郑州 450046)

时滞微分方程在物理学、化学、环境科学和生态学等领域中有着广泛的应用, 但时滞微分方程解析解很难获得, 并且时滞项的存在使得时滞微分方程的理论分析具有一定的困难, 因此研究时滞微分方程的数值解显得十分必要.数值方法的稳定性在时滞微分方程数值解的研究中占有重要地位.在稳定性理论中, 时滞微分方程数值方法的稳定性根据时滞的大小是否对其产生影响, 分为时滞相关稳定性和时滞无关稳定性.时滞微分方程的时滞相关稳定性[1-3], 也称为D-稳定性.数值方法的D-稳定性要求: 数值方法应用于渐近稳定的时滞微分方程

得到的差分方程对所有的步长h=τ/m都要满足渐近稳定的条件.

为了弱化D-稳定性的条件, Hu等[4]提出一种稳定条件比较宽松的时滞相关稳定性-弱时滞相关稳定性, 并且得到了中立型时滞微分方程Runge-Kutta 方法和线性多步法的弱时滞相关稳定性的稳定性判据.弱时滞相关稳定性只需要存在一个正整数m, 使得数值方法应用于一个渐近稳定的时滞微分方程时, 在步长h=τ/m下所得的差分系统满足渐近稳定性条件即可.

本工作研究中立型时滞微分方程

式中: 谱半径ρ(N)<1;y(t)∈Rd; 参数矩阵L,M,N ∈Rd×d; 时滞常数τ >0.

方程(2)数值方法的稳定性研究目前已有很多成果[5-11].Rosenbrock 方法是一类常用于求解刚性常微分方程的数值方法, 该方法计算便捷易行, 具有良好的稳定性.曹学年等[12]最早将Rosenbrock 方法用于求解时滞微分方程, 证明了这类方法是GP-稳定的; 丛玉豪等[13]研究了多维时滞微分方程Rosenbrock 方法的GPL-稳定性; Zhao等[14]研究了中立型时滞微分代数方程Rosenbrock 方法的GP-稳定性; 覃婷婷等[15]构造了含分布时滞的中立型微分方程的Rosenbrock 方法, 获得了一些稳定性的结果; 王艳沛[16]研究了含分布时滞的时滞微分系统Rosenbrock 方法的弱时滞相关稳定性, 并给出了稳定性判据.

本工作在方程(2)渐近稳定的前提下讨论Rosenbrock 方法的弱时滞相关稳定性, 并根据幅角原理获得稳定性判据, 给出的数值例子验证了所得结论的可行性和有效性.

1 理论分析

方程(2)的特征方程为

式中:I为单位矩阵.式(3)的根称为特征根.方程(2)的渐近稳定性由特征根的分布决定, 也就是说, 方程(2)是渐近稳定的当且仅当所有特征根都位于复平面的左半平面内, 即特征根的实部小于0.假如存在不稳定的特征根, 即特征根的实部大于或者等于0, 由下面引理可知, 这些不稳定的特征根在一个有界区域内.

引理1[4]假设ρ(N)<1,s为式(3)中一个不稳定的特征根, 即Res≥0, 则存在一个正定矩阵H, 满足H −N∗HN=I, 且有

式中:

定义1[4]假设引理1 的条件满足, 且边界β已经获得, 那么s平面上Dβ是一个闭集,

Dβ的边界记为Γβ,Γβ包括两部分:

引理2[4]假设引理1 的条件满足, 且边界半径为β, 那么方程(2)是渐近稳定的, 当且仅当

∆ΓβargP(z)表示P(z)沿着Γβ的正方向转动一周后幅角的改变量.

定义2[4]假设给定矩阵L,M,N和时滞τ >0,ρ(N)<1, 方程(2)是渐近稳定的.如果存在正整数m, 使得h=τ/m, 且数值解yn满足

则称数值方法应用于方程(2)时是弱时滞相关稳定的.

2 Rosenbrock 方法的弱时滞相关稳定性

对于非自治系统

应用s级Rosenbrock 方法[13],

式中:Ki是级值的近似;h是步长;tn=t0+nh; J 表示Jacobi 矩阵

对于中立型时滞微分方程

式中:τ >0;y ∈Rd; 函数f为[t0,T]×Rd×Rd×Rd −→Rd.

对方程(6)应用s级Rosenbrock 方法, 可得

式中:h=τ/m;tn=t0+nh;yn是y(tn)的近似值;yn,j,yn−m,Kn,j,Kn−m,j,Kn−m分别是y(tn+cjh),y(tn −τ), ˙y(tn+cjh), ˙y(tn+cjh −τ), ˙y(tn −τ)的近似值;αij=0(i≤j);βij=0(i≤j);γij=0(i

将s级Rosenbrock 方法应用于方程(2), 可得到差分方程

引理3差分方程(8)的特征多项式PRB(z)可以表示为

式中:P=(αij)s×s;Q=(βij)s×s;R=(γij)s×s;⊗表示Kronecker 积;b表示(b1,b2,··· ,bs)T;e表示(1,1,··· ,1)T;Is,Id,Isd分别表示s,d,sd阶单位矩阵.

证明 差分方程(8)可以表示为

式中:

因此, 差分方程(8)的特征多项式PRB(z)可以表示为式(9).

定理1对于s级Rosenbrock 方法, 如果存在正整数m, 使得如下条件成立:

(1)h=τ/m;

(2) 给定矩阵L,M,N和τ,ρ(N)<1, 系统(2)是渐近稳定的, 即引理2 成立;

(3) 矩阵h(R ⊗L)+h(P ⊗L)的所有特征根不等于1;

(4) 特征多项式PRB(z)在复平面的单位圆µ={z:|z|= 1}上不等于0, 且在单位圆上的幅角的变化满足

则s级Rosenbrock 方法应用于方程(2)时是弱时滞相关稳定的.

证明 特征方程PRB(z)=0 中最高次项zm+1的系数矩阵为

在条件(3)满足的情况下, 矩阵Isd −h(R ⊗L)−h(P ⊗L)是非奇异的, 因此特征多项式的最高次项zm+1的系数矩阵是非奇异的.于是, 特征方程PRB(z)=0 有d(s+1)(m+1)个特征根, 如果有重数, 按重数计算.

差分方程(8)渐近稳定的充要条件是特征方程PRB(z) = 0 的所有特征根都在单位圆内.由幅角原理可知, 条件(4)确保了特征方程PRB(z) = 0 满足这一性质, 因此差分方程(8)是渐近稳定的, 即此时的Rosenbrock 方法在弱稳定意义下是时滞相关稳定的.定理得证.

3 数值算例

在如下数值例子中, 选用3 级Rosenbrock 方法, 系数矩阵如下:

例1 考虑2 维中立型时滞微分方程, 其对应的参数矩阵为

根据引理1,ρ(N)=0.85<1, 不稳定域的半径β=18.129 6.

(1) 当τ=1 时, 经过计算, ∆ΓβargP(z)=0, 引理2 的条件满足, 即定理1 中条件(2)成立,此时的时滞微分方程(2)是渐近稳定的.下面验证此时Rosenbrock 方法的稳定性.

(a) 取m= 10,h=τ/m= 0.1, 通过计算, 矩阵h(R ⊗L)+h(P ⊗L)的所有特征根为−0.031 7,−0.031 7,−0.031 7,−0.021 1,−0.021 1,−0.021 1, 所有特征值均不等于1, 即定理1 中条件(3)满足.

图1 τ =1, m=10 时算例1 的数值解Fig.1 Numerical solution of example 1 for τ =1, m=10

(b) 取m= 80,h=τ/m= 0.012 5, 通过计算, 矩阵h(R ⊗L)+h(P ⊗L)的所有特征根为−0.004 0,−0.004 0,−0.004 0,−0.002 6,−0.002 6,−0.002 6, 所有特征值均不等于1, 即定理1 中条件(3)满足.

图2 τ =1, m=80 时算例1 的数值解Fig.2 Numerical solution of example 1 for τ =1, m=80

(2) 当τ= 2.5 时, 经过计算, ∆ΓβargP(z) = 0, 引理2 的条件满足, 即定理1 中条件(2)成立, 则此时的时滞微分系统(2)是渐近稳定的.下面验证数值方法的稳定性.

(a) 取m= 5,h=τ/m= 0.5, 通过计算, 矩阵h(R ⊗L) +h(P ⊗L)的所有特征根为−0.158 5,−0.158 5,−0.158 5,−0.105 7,−0.105 7,−0.105 7, 所有特征值均不等于1, 即定理1 中条件(3)满足.

图3 τ =2.5, m=5 时算例1 的数值解Fig.3 Numerical solution of example 1 for τ =2.5, m=5

(b) 取m= 50,h=τ/m= 0.05, 通过计算, 矩阵h(R ⊗L)+h(P ⊗L)的所有特征根为−0.015 8,−0.015 8,−0.015 8,−0.010 6,−0.010 6,−0.010 6, 所有特征值均不等于1, 即定理1 中条件(3)满足.

图4 τ =2.5, m=50 时算例1 的数值解Fig.4 Numerical solution of example 1 for τ =2.5, m=50

例2 考虑4 维中立型时滞微分方程, 同样选取3 级Rosenbrock 方法, 以2 步2 级Runge-Kutta 方法作对比.微分方程对应的参数矩阵如下:

根据引理1,ρ(N)=0.085 4<1, 不稳定域的半径β=23.378 2.

(1) 当τ= 0.1 时, 经过计算, ∆ΓβargP(z) = 0, 引理2 的条件满足, 即定理1 中条件(2)成立, 此时的时滞微分系统(2)是渐近稳定的.下面验证数值方法的稳定性.

(a) 取m= 5,h=τ/m= 0.02, 通过计算, 矩阵h(R ⊗L)+h(P ⊗L)的所有特征根均不等于1, 即定理1 中条件(3)满足.

图5 τ =0.1, m=5 时算例2 的数值解Fig.5 Numerical solution of example 2 for τ =0.1, m=5

(b) 取m=100,h=τ/m=0.001, 通过计算, 矩阵h(R ⊗L)+h(P ⊗L)的所有特征根均不等于1, 即定理1 中条件(3)满足.

图6 τ =0.1, m=100 时算例2 的数值解Fig.6 Numerical solution of example 2 for τ =0.1, m=100

(2) 当τ= 0.3 时, 经过计算, ∆ΓβargP(z) = 2, 由引理2, 则此时的时滞微分系统(2)不是渐近稳定的.

(a)取m=80,h=τ/m=0.003 75,∆µPRB(z)=1 293.982≈1 294,d(s+1)(m+1)=4(3+1)(80+1)=1 296, 即∆µPRB(z)̸=d(s+1)(m+1), 即定理1 中条件(4)不满足, 那么由定理1 可知, 3 级Rosenbrock 方法应用于该系统是不稳定的, 数值解如图7(a)所示.可以看出, 数值解发散.图7(b)为2 步2 级Runge-Kutta 方法得到的结果, 同样数值解发散.

图7 τ =0.3, m=80 时算例2 的数值解Fig.7 Numerical solution of example 2 for τ =0.3, m=80

4 结束语

本工作研究了Rosenbrock 方法应用于中立型时滞微分方程时的弱时滞相关稳定性, 给出了稳定性的充分条件, 即对于渐近稳定的中立型时滞微分系统, 如果存在合适的正整数m使得定理1 成立, 那么可以通过Rosenbrock 方法得到渐近稳定的数值解.