姚潔怡， 王 琦

(廣東工業(yè)大學 應用數(shù)學學院， 廣東 廣州 510520)

加拿大學者Michael C. Mackey[1]在20世紀70年代對生理系統(tǒng)的研究做了開創(chuàng)性的工作。此后，因為Mackey-Glass系統(tǒng)作為一個基礎的生物學例子，對生理調控也有一定的指導意義，其動力學性質受到越來越多的關注[2-4]。許多學者將他們的研究集中在Mackey-Glass系統(tǒng)的分支、混沌等動力學行為上。實際上，對于具有任何參數(shù)的時滯非線性動力系統(tǒng)，當時滯增加超過某個臨界值時都可能導致系統(tǒng)從穩(wěn)定平衡變得不穩(wěn)定，并導致分支。不同時滯下的穩(wěn)定域如何變化？分支發(fā)生時如何確定時滯的臨界值？質量行為如何依賴于時滯的形式？這些問題促使我們研究時滯變化對系統(tǒng)動力學行為的影響。考慮到科學計算和實時仿真的需要，我們的興趣集中在原系統(tǒng)對應的離散動力學系統(tǒng)的行為上[5-9]。在大多數(shù)情況下，希望由微分方程導出的差分方程保持相應連續(xù)時間模型的動態(tài)特征。即離散時間模型與連續(xù)時間模型“動態(tài)一致”。Ford和Wulf[10]運用根軌跡法，證明了線性多步法的收斂階與原時滯微分方程相同，其穩(wěn)定區(qū)域收斂于原方程的穩(wěn)定區(qū)域；2000年，他們又用歐拉法研究了一類帶參數(shù)的時滯微分方程并且證明了離散格式與連續(xù)時間模型“動態(tài)一致”[11]。這意味著，對于所有足夠小的步長，歐拉法離散模型經(jīng)歷與相應的連續(xù)時間模型相同類型的Hopf分支。

近年來，研究人員發(fā)現(xiàn)許多差分方法都是能夠保持動態(tài)一致的。利用時滯τ作為參數(shù)，DING Xiao-hua等[12]應用梯形方法研究了Mackey-Glass系統(tǒng)的動力性。SU Huan等[13-14]用非標準有限差分方法研究了Mackey-Glass系統(tǒng)的Hopf分支的存在性與穩(wěn)定性并且計算了分支方向，對Hopf分支的控制作了分析。然而，用改進歐拉法求解時滯微分方程的文獻還很少，改進歐拉法能否保持時滯微分方程的動力性有待于研究。在本文中，利用改進歐拉法對Mackey-Glass系統(tǒng)進行離散，研究離散化后的系統(tǒng)在正平衡點處的穩(wěn)定情況和Hopf分支存在性。

1 正平衡點的穩(wěn)定性與局部Hopf分支

對于時滯微分方程

(1)

其中p(t)表示血液循環(huán)系統(tǒng)中t時刻成熟細胞的密度，τ是時間延遲，t≥0，β、θ、n、γ都是正常數(shù)。

令p(t)=θx(t)，方程(1)變成

(2)

這里a=β/θ，讓u(t)=x(τt)，則方程(2)可以寫成

(3)

下面將改進歐拉法應用于方程(3)，取步長h=1/m，令

建立預報-校正系統(tǒng)

(4)

其中un=φ(nh)，-m≤n≤0，且

(5)

un為φ(nh)的近似值。用un-m來近似時滯項u(tn-τ)。

將式(5)代入式(4)，有

化簡得到改進歐拉法遞推格式

(6)

假設u*是方程的一個不動點，u*滿足以下方程

γun+1+γu-a=0，

令

F(x)=γxn+1+γx-a，

對所有x≥0，F(xiàn)′(x)=(n+1)γxn+γ>0，F(xiàn)(0)=-a<0，故方程(1)有唯一正不動點u*。

讓yn=un-u*，那么，有

(7)

引進一個新的變量

Yn=(yn,yn-1,…,yn-m)T，

那么存在一個映射

Yn+1=F(Yn,τ)，

其中

F=(F0,F1,…,Fm)T，

(8)

顯然，原點是上式的一個不動點，并且它的線性部分為

Yk+1=AYk，

因此，A的特征方程為

(9)

引理1 對于充分小的正數(shù)τ，方程(9)的所有根都小于1。

證明當τ=0時，方程(9)等價于λm+1-λm=0，方程有一個m重根和一個單根λ=1。

考慮方程(9)的根λ(τ)，使得λ(0)=1，方程(9)是關于τ可微的，則有

(10)

(11)

所以，隨著τ>0，沒有特征根穿過λ=1。因此，對于充分小的正數(shù)τ，方程(9)的所有特征根都在單位圓內。

隨著τ的變化，當方程(9)出現(xiàn)一對共軛的特征根穿過單位圓時，Hopf分支產生。要證明Hopf分支的存在性，需要找到單位圓上ω的值。

用eiω表示單位圓上的根，則當ω∈(0,π]時，eiω是方程(9)的根，當且僅當

(12)

分離實部和虛部，有

(13)

(14)

所以

(15)

(16)

其中[·]為取整函數(shù)。

(17)

其中ωi滿足方程(16)。

證明根據(jù)方程(14)和(15)，有

得證。

由引理4可得方程(9)零解的穩(wěn)定性。所以關于方程(6)我們有如下定理。

2 數(shù)值實驗

圖1 步長h=1/2時方程(3)的數(shù)值解

圖2 步長h=1/4時方程(3)的數(shù)值解

圖3 步長h=1/8時方程(3)的數(shù)值解

根據(jù)定理1和表1，由圖1—圖4可見，當τ∈[0,τ0)時，方程(3)的數(shù)值解都是漸近穩(wěn)定的；而當τ>τ0時，數(shù)值解不穩(wěn)定。數(shù)值現(xiàn)象與理論分析結果是一致的。

表1 不同步長下分支點的值

圖4 步長h=1/16時方程(3)的數(shù)值解

3 小結

通過與文獻[16]中給出的原系統(tǒng)狀態(tài)圖與動力學性質比較，改進歐拉法的數(shù)值模擬能夠正確地反應原系統(tǒng)的動力學性質，即保持其動態(tài)一致性。