張明坤 王艷沛 趙歡歡

(1.上海大學(xué)理學(xué)院, 上海 200444;2.鄭州航空工業(yè)管理學(xué)院數(shù)學(xué)學(xué)院, 鄭州 450046)

時滯微分方程在物理學(xué)、化學(xué)、環(huán)境科學(xué)和生態(tài)學(xué)等領(lǐng)域中有著廣泛的應(yīng)用, 但時滯微分方程解析解很難獲得, 并且時滯項的存在使得時滯微分方程的理論分析具有一定的困難, 因此研究時滯微分方程的數(shù)值解顯得十分必要.數(shù)值方法的穩(wěn)定性在時滯微分方程數(shù)值解的研究中占有重要地位.在穩(wěn)定性理論中, 時滯微分方程數(shù)值方法的穩(wěn)定性根據(jù)時滯的大小是否對其產(chǎn)生影響, 分為時滯相關(guān)穩(wěn)定性和時滯無關(guān)穩(wěn)定性.時滯微分方程的時滯相關(guān)穩(wěn)定性[1-3], 也稱為D-穩(wěn)定性.數(shù)值方法的D-穩(wěn)定性要求: 數(shù)值方法應(yīng)用于漸近穩(wěn)定的時滯微分方程

得到的差分方程對所有的步長h=τ/m都要滿足漸近穩(wěn)定的條件.

為了弱化D-穩(wěn)定性的條件, Hu等[4]提出一種穩(wěn)定條件比較寬松的時滯相關(guān)穩(wěn)定性-弱時滯相關(guān)穩(wěn)定性, 并且得到了中立型時滯微分方程Runge-Kutta 方法和線性多步法的弱時滯相關(guān)穩(wěn)定性的穩(wěn)定性判據(jù).弱時滯相關(guān)穩(wěn)定性只需要存在一個正整數(shù)m, 使得數(shù)值方法應(yīng)用于一個漸近穩(wěn)定的時滯微分方程時, 在步長h=τ/m下所得的差分系統(tǒng)滿足漸近穩(wěn)定性條件即可.

本工作研究中立型時滯微分方程

式中: 譜半徑ρ(N)<1;y(t)∈Rd; 參數(shù)矩陣L,M,N ∈Rd×d; 時滯常數(shù)τ >0.

方程(2)數(shù)值方法的穩(wěn)定性研究目前已有很多成果[5-11].Rosenbrock 方法是一類常用于求解剛性常微分方程的數(shù)值方法, 該方法計算便捷易行, 具有良好的穩(wěn)定性.曹學(xué)年等[12]最早將Rosenbrock 方法用于求解時滯微分方程, 證明了這類方法是GP-穩(wěn)定的; 叢玉豪等[13]研究了多維時滯微分方程Rosenbrock 方法的GPL-穩(wěn)定性; Zhao等[14]研究了中立型時滯微分代數(shù)方程Rosenbrock 方法的GP-穩(wěn)定性; 覃婷婷等[15]構(gòu)造了含分布時滯的中立型微分方程的Rosenbrock 方法, 獲得了一些穩(wěn)定性的結(jié)果; 王艷沛[16]研究了含分布時滯的時滯微分系統(tǒng)Rosenbrock 方法的弱時滯相關(guān)穩(wěn)定性, 并給出了穩(wěn)定性判據(jù).

本工作在方程(2)漸近穩(wěn)定的前提下討論Rosenbrock 方法的弱時滯相關(guān)穩(wěn)定性, 并根據(jù)幅角原理獲得穩(wěn)定性判據(jù), 給出的數(shù)值例子驗證了所得結(jié)論的可行性和有效性.

1 理論分析

方程(2)的特征方程為

式中:I為單位矩陣.式(3)的根稱為特征根.方程(2)的漸近穩(wěn)定性由特征根的分布決定, 也就是說, 方程(2)是漸近穩(wěn)定的當(dāng)且僅當(dāng)所有特征根都位于復(fù)平面的左半平面內(nèi), 即特征根的實部小于0.假如存在不穩(wěn)定的特征根, 即特征根的實部大于或者等于0, 由下面引理可知, 這些不穩(wěn)定的特征根在一個有界區(qū)域內(nèi).

引理1[4]假設(shè)ρ(N)<1,s為式(3)中一個不穩(wěn)定的特征根, 即Res≥0, 則存在一個正定矩陣H, 滿足H ?N?HN=I, 且有

式中:

定義1[4]假設(shè)引理1 的條件滿足, 且邊界β已經(jīng)獲得, 那么s平面上Dβ是一個閉集,

Dβ的邊界記為Γβ,Γβ包括兩部分:

引理2[4]假設(shè)引理1 的條件滿足, 且邊界半徑為β, 那么方程(2)是漸近穩(wěn)定的, 當(dāng)且僅當(dāng)

?ΓβargP(z)表示P(z)沿著Γβ的正方向轉(zhuǎn)動一周后幅角的改變量.

定義2[4]假設(shè)給定矩陣L,M,N和時滯τ >0,ρ(N)<1, 方程(2)是漸近穩(wěn)定的.如果存在正整數(shù)m, 使得h=τ/m, 且數(shù)值解yn滿足

則稱數(shù)值方法應(yīng)用于方程(2)時是弱時滯相關(guān)穩(wěn)定的.

2 Rosenbrock 方法的弱時滯相關(guān)穩(wěn)定性

對于非自治系統(tǒng)

應(yīng)用s級Rosenbrock 方法[13],

式中:Ki是級值的近似;h是步長;tn=t0+nh; J 表示Jacobi 矩陣

對于中立型時滯微分方程

式中:τ >0;y ∈Rd; 函數(shù)f為[t0,T]×Rd×Rd×Rd ?→Rd.

對方程(6)應(yīng)用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 方法應(yīng)用于方程(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 方法, 如果存在正整數(shù)m, 使得如下條件成立:

(1)h=τ/m;

(2) 給定矩陣L,M,N和τ,ρ(N)<1, 系統(tǒng)(2)是漸近穩(wěn)定的, 即引理2 成立;

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

(4) 特征多項式PRB(z)在復(fù)平面的單位圓μ={z:|z|= 1}上不等于0, 且在單位圓上的幅角的變化滿足

則s級Rosenbrock 方法應(yīng)用于方程(2)時是弱時滯相關(guān)穩(wěn)定的.

證明 特征方程PRB(z)=0 中最高次項zm+1的系數(shù)矩陣為

在條件(3)滿足的情況下, 矩陣Isd ?h(R ?L)?h(P ?L)是非奇異的, 因此特征多項式的最高次項zm+1的系數(shù)矩陣是非奇異的.于是, 特征方程PRB(z)=0 有d(s+1)(m+1)個特征根, 如果有重數(shù), 按重數(shù)計算.

差分方程(8)漸近穩(wěn)定的充要條件是特征方程PRB(z) = 0 的所有特征根都在單位圓內(nèi).由幅角原理可知, 條件(4)確保了特征方程PRB(z) = 0 滿足這一性質(zhì), 因此差分方程(8)是漸近穩(wěn)定的, 即此時的Rosenbrock 方法在弱穩(wěn)定意義下是時滯相關(guān)穩(wěn)定的.定理得證.

3 數(shù)值算例

在如下數(shù)值例子中, 選用3 級Rosenbrock 方法, 系數(shù)矩陣如下:

例1 考慮2 維中立型時滯微分方程, 其對應(yīng)的參數(shù)矩陣為

根據(jù)引理1,ρ(N)=0.85<1, 不穩(wěn)定域的半徑β=18.129 6.

(1) 當(dāng)τ=1 時, 經(jīng)過計算, ?ΓβargP(z)=0, 引理2 的條件滿足, 即定理1 中條件(2)成立,此時的時滯微分方程(2)是漸近穩(wěn)定的.下面驗證此時Rosenbrock 方法的穩(wěn)定性.

(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 的數(shù)值解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 的數(shù)值解Fig.2 Numerical solution of example 1 for τ =1, m=80

(2) 當(dāng)τ= 2.5 時, 經(jīng)過計算, ?ΓβargP(z) = 0, 引理2 的條件滿足, 即定理1 中條件(2)成立, 則此時的時滯微分系統(tǒng)(2)是漸近穩(wěn)定的.下面驗證數(shù)值方法的穩(wěn)定性.

(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 的數(shù)值解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 的數(shù)值解Fig.4 Numerical solution of example 1 for τ =2.5, m=50

例2 考慮4 維中立型時滯微分方程, 同樣選取3 級Rosenbrock 方法, 以2 步2 級Runge-Kutta 方法作對比.微分方程對應(yīng)的參數(shù)矩陣如下:

根據(jù)引理1,ρ(N)=0.085 4<1, 不穩(wěn)定域的半徑β=23.378 2.

(1) 當(dāng)τ= 0.1 時, 經(jīng)過計算, ?ΓβargP(z) = 0, 引理2 的條件滿足, 即定理1 中條件(2)成立, 此時的時滯微分系統(tǒng)(2)是漸近穩(wěn)定的.下面驗證數(shù)值方法的穩(wěn)定性.

(a) 取m= 5,h=τ/m= 0.02, 通過計算, 矩陣h(R ?L)+h(P ?L)的所有特征根均不等于1, 即定理1 中條件(3)滿足.

圖5 τ =0.1, m=5 時算例2 的數(shù)值解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 的數(shù)值解Fig.6 Numerical solution of example 2 for τ =0.1, m=100

(2) 當(dāng)τ= 0.3 時, 經(jīng)過計算, ?ΓβargP(z) = 2, 由引理2, 則此時的時滯微分系統(tǒng)(2)不是漸近穩(wěn)定的.

(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 方法應(yīng)用于該系統(tǒng)是不穩(wěn)定的, 數(shù)值解如圖7(a)所示.可以看出, 數(shù)值解發(fā)散.圖7(b)為2 步2 級Runge-Kutta 方法得到的結(jié)果, 同樣數(shù)值解發(fā)散.

圖7 τ =0.3, m=80 時算例2 的數(shù)值解Fig.7 Numerical solution of example 2 for τ =0.3, m=80

4 結(jié)束語

本工作研究了Rosenbrock 方法應(yīng)用于中立型時滯微分方程時的弱時滯相關(guān)穩(wěn)定性, 給出了穩(wěn)定性的充分條件, 即對于漸近穩(wěn)定的中立型時滯微分系統(tǒng), 如果存在合適的正整數(shù)m使得定理1 成立, 那么可以通過Rosenbrock 方法得到漸近穩(wěn)定的數(shù)值解.