姚潔怡， 王 琦

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

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

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

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

對于時(shí)滯微分方程

(1)

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

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

(2)

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

(3)

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

建立預(yù)報(bào)-校正系統(tǒng)

(4)

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

(5)

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

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

化簡得到改進(jìn)歐拉法遞推格式

(6)

假設(shè)u*是方程的一個(gè)不動(dòng)點(diǎn)，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)有唯一正不動(dòng)點(diǎn)u*。

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

(7)

引進(jìn)一個(gè)新的變量

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

那么存在一個(gè)映射

Yn+1=F(Yn,τ)，

其中

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

(8)

顯然，原點(diǎn)是上式的一個(gè)不動(dòng)點(diǎn)，并且它的線性部分為

Yk+1=AYk，

因此，A的特征方程為

(9)

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

證明當(dāng)τ=0時(shí)，方程(9)等價(jià)于λm+1-λm=0，方程有一個(gè)m重根和一個(gè)單根λ=1。

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

(10)

(11)

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

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

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

(12)

分離實(shí)部和虛部，有

(13)

(14)

所以

(15)

(16)

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

(17)

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

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

得證。

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

2 數(shù)值實(shí)驗(yàn)

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

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

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

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

表1 不同步長下分支點(diǎn)的值

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

3 小結(jié)

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