具分段常數(shù)變量時(shí)滯造血模型的分支分析

陳斯養(yǎng)，黃曉宇

（陜西師范大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，西安710062）

帶有分段常數(shù)變量的微分方程是近年來引起廣泛關(guān)注的一類泛函微分方程.這類方程是連續(xù)和離散變量的混合體，具有微分方程和差分方程的雙重性質(zhì)，在生態(tài)數(shù)學(xué)中具有重要的應(yīng)用.文獻(xiàn)［1］研究了具分段常數(shù)變量微分方程的全局吸引性，文獻(xiàn)［2－4］討論了具分段常數(shù)變量微分方程的振動及非振動性，文獻(xiàn)［5－6］利用概周期型序列研究了具有分段常數(shù)變量微分方程的概周期解、漸近概周期解和偽概周期解的存在性和唯一性，文獻(xiàn)［7］用中心流形理論分析了具有分段常數(shù)變量的捕食－被捕食模型的Flip分支、Hopf分支及Neimark－Sacker分支.

本文給出了具分段常數(shù)變量的時(shí)滯造血模型其中，P（t）表示血液循環(huán)過程中成熟干細(xì)胞的密度，l，m是骨髓中未成熟干細(xì)胞由產(chǎn)生到成為成熟干細(xì)胞所需時(shí)間，β，b，α，r均為正參數(shù).

為了便于對模型（1）進(jìn)行研究，在其等式兩邊同乘以ert，并分別在區(qū)間［n，t）上積分，整理可得

其中，t∈［n，n＋1）（n＝0，1，2，…），令t→n＋1，可得時(shí)滯差分方程模型

本文討論模型（1）的等價(jià)模型（3）當(dāng)1≤m≤l≤2時(shí)正平衡態(tài)的局部漸近穩(wěn)定性，并對其Flip分支和Neimark－Sacker（以下均簡寫為N－S）分支進(jìn)行分支分析.

1 正平衡態(tài)的存在性及局部穩(wěn)定性

為了分析（3）正平衡態(tài)?P的局部穩(wěn)定性，令X（n）＝P（n）－，則（3）等價(jià)轉(zhuǎn)化為：

將（4）在X＝0處進(jìn)行Taylor展開，可得其線性近似系統(tǒng)為：

由（5）得特征方程為：

由差分方程的穩(wěn)定性理論可知，當(dāng)且僅當(dāng)所有的特征根滿足｜λi｜＜1，i＝1，2，3，…，l＋1，正平衡態(tài)是局部漸近穩(wěn)定的.

根據(jù)Jury判據(jù)［8］，可以得到以下關(guān)于正平衡態(tài)?P局部漸近穩(wěn)定的充要條件.

定理1對模型（3），以下結(jié)論成立：

（Ⅰ）當(dāng)m＝l＝1時(shí)，?P局部漸近穩(wěn)定的充要條件為：

（Ⅱ）當(dāng)m＝l＝2時(shí)，?P局部漸近穩(wěn)定的充要條件為：

（Ⅲ）當(dāng)m＝1，l＝2時(shí)，

2 分支的存在性

由分支理論［9］可知，對應(yīng)于特征值λ＝－1的分支稱為Flip分支，對應(yīng)于λ1，2＝e±iθ，0＜θ＜π的分支稱為N－S分支.這一節(jié)，選取β0作為分支參數(shù)討論（3）的Flip分支和N－S分支的存在性條件.

對于特征方程（6），記G（λ）＝λl＋1－e－rλl－（1

證明若定理?xiàng)l件成立，則G（λ）＝λ2－e－rλ＋1，此時(shí)，對于方程G（λ）＝0，求解可得其存在一對共軛復(fù)根且滿足｜λ1，2｜＝1.故（3）在產(chǎn)生N－S分支.

證明若定理?xiàng)l件成立，則G（λ）＝λ3－e－rλ2－（）此時(shí)，方程Gλ＝0可分解為

且滿足｜λ3｜＜1.故（3）在產(chǎn)生N－S分支.

證明若定理的條件①成立，則G（λ）＝λ3－e－rλ2－（1－e－r）λ＋2e－r，此時(shí)，方程G（λ）＝0因式分解為G（λ）＝（λ＋1）（λ2－（1＋e－r）λ＋2e－r），可得λ1＝－1，｜λ2，3｜＜1，故（3）在產(chǎn)生Flip分支；

若定理的條件②成立，則G（λ）＝λ3－e－rλ2－，此時(shí)，對方程G（λ）＝0作因式分解

注1當(dāng)定理2、定理3的條件及定理4的條件②成立時(shí)，模型（3）不存在Flip分支；當(dāng)定理4的條件①成立時(shí)，模型（3）不存在N－S分支.

3 分支的方向和穩(wěn)定性

這一節(jié)利用文獻(xiàn)［9］中的規(guī)范化理論和中心流形的計(jì)算方法研究和討論了分支的方向及其穩(wěn)定性.現(xiàn)將模型（3）等價(jià)轉(zhuǎn)化為l＋1維差分方程組

為了研究（3）的分支方向及穩(wěn)定性，需將（7）轉(zhuǎn)化為ρ（n＋1）＝Uρ（n）＋Q（ρ（n））.其中，ρ（n）＝（P1，P2，…，Pm＋1，…，Pl＋1）∈Rl＋1，U為l＋1維方陣，Q（ρ）＝O（‖ρ‖2）是光滑函數(shù)，在ρ＝0處，Q（ρ）的Taylor展開為，其中E和F是多重線性函數(shù)，對平面向量x＝（x1，x2，…，xl＋1）T，y＝（y1，y2，…，yl＋1）T和z＝（z1，z2，…，zl＋1）T取值，在坐標(biāo)下的分量分別為

其中，i＝1，2，3，…，l＋1.

3.1 當(dāng)m＝l＝1，2時(shí)分支的方向和穩(wěn)定性

當(dāng)m＝l＝2時(shí)，（3）轉(zhuǎn)化為等價(jià)差分方程組若定理3的條件成立，則模型（8）在正平衡態(tài)（，）T產(chǎn)生N－S分支.

作變量變換P1（n）＝＋ξ1（n），P2（n）＝＋ξ2（n），P3（n）＝＋ξ3（n），將（8）轉(zhuǎn)化為

將（9）的第一個(gè)方程在ξ2（n）＝0進(jìn)行Taylor展開，得到等價(jià)方程組

此系統(tǒng)可表示為：

其中

系統(tǒng)（10）在中心流形上的限制規(guī)范化方程為z（n＋1）＝eiθz（n）（1＋l｜z（n）｜2）＋O（｜z（n）｜4），其中，z（n）＝〈p1，ξ（n）〉，而l0＝Rel為N－S分支的臨界規(guī)范形系數(shù)，決定了分支的方向及穩(wěn)定性，且

當(dāng)m＝l＝1時(shí)，若定理2的條件成立，則運(yùn)用同樣的方法可得N－S分支的臨界規(guī)范形系數(shù)為：

定理5若定理2的條件成立，當(dāng)l0、l1＜0（＞0）時(shí)，模型（3）存在唯一穩(wěn)定（不穩(wěn)定）的超（亞）臨界的N－S分支.

3.2 當(dāng)m＝1，l＝2時(shí)分支的方向和穩(wěn)定性

當(dāng)m＝1，l＝2時(shí)，（3）轉(zhuǎn)化為等價(jià)差分方程組

將（12）的第一個(gè)方程在（ξ2（n），ξ3（n））＝（0，0）進(jìn)行Taylor展開，得到等價(jià)方程組

此系統(tǒng)可表示為

其中，

而

其中，B1＝（α－1）bα－（α＋1）（）α，B2＝（α2－3α＋2）b2α－4（α2－1）bα（）α＋（α2＋3α＋2）（）2α.

向量q2～（－1，1，－1）T，p2～分別為矩陣B0和的特征向量，即滿足B0q2＝e，p2＝e－iθp2.為了使向量q2和p2滿足〈p2，q2〉＝p21q21＋p22q22＋p23q23＝1，可取q2＝（－1，1，－1）T，p2＝.多重線性函數(shù)E和F分別取值為

而

系統(tǒng)（13）在中心流形上的限制規(guī)范的化方程為ω（n＋1）＝－ω（n）＋d0ω3（n）＋O（ω4（n）），其中，ω（n）＝〈p2，ξ（n）〉，而d0為Flip分支的臨界規(guī)范形系數(shù)，決定了Flip分支的方向及穩(wěn)定性.計(jì)算得該系數(shù)為：

其中，

定理6若定理4的條件①成立，當(dāng)d0＜0（＞0）時(shí)，模型（9）存在唯一穩(wěn)定（不穩(wěn)定）的、超（亞）臨界的Flip分支.

若定理4的條件②成立，則模型（11）在正平衡態(tài)（?P，?P，?P）T產(chǎn)生N－S分支.作變量變換以及利用Taylor展開公式，可得到系統(tǒng)（11）的等價(jià)方程組

此系統(tǒng)可表示為：

其中，

而

多重線性函數(shù)E和F分別取值為：

可計(jì)算得N－S分支的臨界規(guī)范形系數(shù)為：

定理7若定理4的條件②成立，當(dāng)l2＜0（＞0）時(shí)，模型（3）存在唯一穩(wěn)定（不穩(wěn)定）的、超（亞）臨界的N－S分支.

4 數(shù)值模擬

這一部分對模型（3）進(jìn)行了數(shù)值模擬，通過圖形對模型（3）進(jìn)行分支分析并體現(xiàn)了該生物系統(tǒng)復(fù)雜的動力學(xué)行為，如周期倍分和混沌現(xiàn)象等.下面通過實(shí)例驗(yàn)證文中所得定理與數(shù)值計(jì)算的一致性.

例1在模型（3）或其等價(jià)方程組（8）中取m＝l＝1，r＝0.5，b＝1，α＝4，根據(jù)定理2，可得分支參數(shù)的臨界值為β0＝4.362，此時(shí)唯一的正平衡態(tài)為P＝1.6671.于是，當(dāng)β＜β0時(shí)，P為局部漸近穩(wěn)定的，當(dāng)β＞β0時(shí)，P變?yōu)椴环€(wěn)定.在臨界參數(shù)值β0，有特征值λ1，2＝e±iθ，0＜θ＜π，此時(shí)，N－S分支產(chǎn)生，見β－P分支圖1.計(jì)算可得分支的臨界規(guī)范形系數(shù)l1＜0，由定理5，可知當(dāng)β＞β0時(shí)，從P分支出唯一穩(wěn)定的N－S分支.

當(dāng)β＝4.2＜β0時(shí)，P為局部漸近穩(wěn)定的，此時(shí)，所有的解均趨于正平衡態(tài)P＝P1＝1.6493（見圖2）；當(dāng)β＝4.66＞β0時(shí)，在?P＝?P2＝1.6984產(chǎn)生N－S分支及分支周期解（見圖3）.

圖1 β－P分支圖Fig.1 Branch map about rand P

圖2 穩(wěn)定圖（β＝4.2＜β0）Fig.2 Stability map（β＝4.2＜β0）

圖3 N－S分支及周期解圖（β＝4.66＞β0）Fig.3 Branch map of N－S and periodic solution（β＝4.66＞β0）

例2當(dāng)m＝l＝2時(shí)，在模型（3）或其等價(jià)方程組（8）中取r＝3，b＝1，α＝3，根據(jù)定理3，可計(jì)算得分支參數(shù)的臨界值為β0＝9.2452，此時(shí)，模型（3）存在唯一的正平衡態(tài)P＝1.2769.當(dāng)β＜β0時(shí)，正平衡態(tài)為局部漸近穩(wěn)定的，當(dāng)β＞β0時(shí)，P變?yōu)椴环€(wěn)定.在臨界參數(shù)值β0，有特征值λ1，2＝e±iθ，0＜θ＜π，此時(shí)，N－S分支產(chǎn)生，見β－P分支圖4.計(jì)算可得分支的臨界規(guī)范形系數(shù)l0＜0，由定理5，可知當(dāng)β＞β0時(shí)，在產(chǎn)生唯一穩(wěn)定的、超臨界的N－S分支.

當(dāng)β＝8.995＜β0時(shí)，＝＝1.2596為局部漸近穩(wěn)定的（見圖5）；當(dāng)β＝9.44＞β0時(shí)，在正平衡態(tài)＝＝1.29產(chǎn)生N－S分支及分支周期解（見圖6）.

圖4 β－P分支圖Fig.4 Branch map about rand P

圖5 穩(wěn)定圖（β＝8.995＜β0）Fig.5 Stability map（β＝8.995＜β0）

圖6 N－S分支及周期解圖（β＝9.44＞β0）Fig.6 Branch map of N－S and periodic solution（β＝9.44＞β0）

例3當(dāng)m＝1，l＝2時(shí)，在模型（3）或其等價(jià)方程組（17）中取r＝0.5，b＝1，α＝3，根據(jù)定理4的條件②，可計(jì)算得分支參數(shù)的臨界值為β0＝2.2281，此時(shí)，模型（9）唯一的正平衡態(tài)為?P＝1.5119，對應(yīng)的β－P分支圖（見圖7）.由圖7可知，當(dāng)β＜β0時(shí)，正平衡態(tài)為局部漸近穩(wěn)定的，當(dāng)β＞β0時(shí)，變?yōu)椴环€(wěn)定，在產(chǎn)生N－S分支.計(jì)算可得分支的臨界規(guī)范形系數(shù)l2＜0，由定理7，可知當(dāng)β＞β0時(shí)，在產(chǎn)生唯一穩(wěn)定的N－S分支.

當(dāng)β＝2.1＜β0時(shí)，＝＝1.4736為局部漸近穩(wěn)定的（見圖8）；當(dāng)β＝2.4＞β0時(shí)，在正平衡態(tài)＝＝1.5605產(chǎn)生N－S分支及分支周期解（見圖9）.

圖7 β－P分支圖Fig.7 Branch map about rand P

圖8 穩(wěn)定圖（β＝2.1＜β0）Fig.8 Stability map（β＝2.1＜β0）

圖9 N－S分支及周期解圖（β＝2.4＞β0）Fig.9 Branch map of N－S and periodic solution（β＝2.4＞β0）

例4當(dāng)m＝1，l＝2時(shí)，在模型（3）或其等價(jià)方程組（11）中取r＝1.5，b＝1，α＝1，根據(jù)定理4的條件①，可計(jì)算得分支參數(shù)的臨界值為β0＝3.5247，唯一的正平衡態(tài)為P＊＝1.3498.于是，當(dāng)β＜β0時(shí)，P＊為局部漸近穩(wěn)定的，當(dāng)β＞β0時(shí)，P＊變?yōu)椴环€(wěn)定.在臨界參數(shù)值β＝β0，有特征值λ＝－1，此時(shí)Flip分支產(chǎn)生.計(jì)算分支的臨界規(guī)范形系數(shù)可得d0＜0，由定理6，可知當(dāng)β＞β0時(shí)，在P＊產(chǎn)生唯一穩(wěn)定的Flip分支.

當(dāng)β＝3.36＜β0時(shí)，P＊為局部漸近穩(wěn)定的，此時(shí)P＊＝P7＊＝1.24（見圖10）；當(dāng)β＝3.68＞β0時(shí)，在P＊＝P8＊＝1.4533產(chǎn)生Flip分支及分支周期解（見圖11）.

圖10 穩(wěn)定圖（β＝3.36＜β0）Fig.10 Stability map（β＝3.36＜β0）

圖11 Flip分支及周期解圖（β＝3.68＞β0）Fig.11 Branch map of Flip and periodic solution（β＝3.68＞β0）

［1］ Wang Youbin，Yan Jurang.Necessary and Sufficient Condition for the Global Attractivity of the Trivial Solution of a Delay equation with Continuous and Piecewise Constant Arguments［J］.Appl Math.Lett，1997，10（5）：91－96.

［2］ Wang Youbin，Yan Jurang.Necessary and Sufficient Conditions for the Oscillation of a Delay Logistic equation with Continuous and Piecewise Constant Arguments［J］.Ann of Diff Eqs，2005，21（3）：435－438.

［3］ Shen J H，Stavroulakis I P.Oscillatory and Nonoscillatory Delay equations with Piecewise Constant Argument［J］.Journal of Mathematical Analysis and Applications，2000，248：385－401.

［4］ 王幼斌，燕居讓.具逐段常變量時(shí)滯微分方程的振動性［J］.科學(xué)通報(bào)，1996，41（5）：393－395.

［5］ AIonso A I，Hong J L，Obay A R.Almost Periodic Type Solutions of Differential equations with Piecewise Constant Argument via Almost Periodic Type Sequences［J］.Appl Math Lett，2000，13：131－137.

［6］ 樸大雄.具逐段常變量的中立型時(shí)滯微分方程的偽概周期解［J］.北京大學(xué)學(xué)報(bào)：自然科學(xué)版，2001，37（3）：298－304.

［7］ 陳斯養(yǎng)，張 艷.具有分段常數(shù)變量的捕食—被捕食模型的分支分析［J］.蘭州大學(xué)學(xué)報(bào)：自然科學(xué)版，2012，48（3）：104－112.

［8］ Mustafa R S K，Merino O.Discrete Dynamical Systems and Difference Equations with Mathematica［M］.New York：A CRC Press Company，2002.

［9］ Yuri A K.Elements of Applied Bifurcation Theory［M］.3rd ed.New York：Springer－Verlag，2004.