用最優(yōu)化方法計(jì)算時(shí)滯微分方程的周期解

馮 蘭，徐 旭

（吉林大學(xué) 數(shù)學(xué)學(xué)院，長春130012）

時(shí)滯微分方程是一類重要的數(shù)學(xué)模型，應(yīng)用廣泛，目前已有許多研究結(jié)果.如時(shí)滯微分方程的數(shù)值穩(wěn)定性分析、時(shí)滯微分方程周期解的存在性、唯一性及其數(shù)值解的分析等［1－5］.時(shí)滯微分方程是復(fù)雜的無窮維系統(tǒng)，系統(tǒng)的平衡點(diǎn)經(jīng)過Hopf分岔可以產(chǎn)生周期解［3］.時(shí)滯系統(tǒng)初值定義在一段區(qū)間上，因此需要尋找合適的初始函數(shù)定位周期解；計(jì)算復(fù)雜性較高，而高效計(jì)算時(shí)滯微分方程周期解的方法目前報(bào)道較少［6－7］.時(shí)滯方程周期解的計(jì)算主要有如下幾種方法：

1）通過Fourier級數(shù)逼近定位周期解［4－5］；

2）通過攝動法近似計(jì)算周期解的近似解析表達(dá)式［8－11］；

3）打靶法［4，12－13］.

Fourier級數(shù)方法的效率較高，但無法確定周期解的穩(wěn)定性；利用攝動法近似計(jì)算周期解一般只針對具有弱非線性項(xiàng)的系統(tǒng)才有效；打靶法雖然可以通過確定Floquet乘子確定周期解的穩(wěn)定性，但工作量大、效率低.目前討論時(shí)滯微分方程周期解的研究結(jié)果較多，但大多數(shù)都是定性分析，主要目的是得到分岔周期解的分岔點(diǎn)、周期及用中心流形和規(guī)范性理論計(jì)算分岔周期解的穩(wěn)定性.在數(shù)值計(jì)算時(shí)滯微分方程周期解的方法中，Newton－Picard法通過降低Jacobi矩陣中非零元素的個(gè)數(shù)提高計(jì)算效率［4－5］.本文方法與Newton－Picard的思路不同，應(yīng)用最優(yōu)化方法求解時(shí)滯系統(tǒng)的周期解，提出了用函數(shù)擬合方法確定初始函數(shù)的方法.結(jié)果表明，與工程上常用的通過函數(shù)插值方法逼近初始函數(shù)的方法相比［4］，本文方法極大減小了計(jì)算量，提高了計(jì)算效率.

1 時(shí)滯方程周期解的計(jì)算

考慮如下具有時(shí)滯τ的微分方程系統(tǒng)：

其中：x是n維向量；τ≥0表示時(shí)滯；φ表示定義在［－τ，0）上的初始函數(shù)；F（·）表示n維實(shí)值向量函數(shù).設(shè)方程（1）的解片段為xt（φ）＝x（t＋θ，φ），這里xt（φ）（θ）∈C（［－τ，0］，?n），其中：－τ≤θ≤0；C是［－τ，0］→?n的連續(xù)函數(shù)Banach空間.

時(shí)滯方程的初始值定義在一個(gè)初始函數(shù)φ（θ）（－τ≤θ≤0）上，系統(tǒng)的周期解應(yīng)滿足xT（φ）＝φ，這里T是周期.為了尋找方程（1）的周期T解，先考慮如下Poincaré映射：

這里xT（φ）為定義在［T－τ，T］上的解片段.于是，計(jì)算時(shí)滯方程的周期解即為尋找映射（2）的固定點(diǎn)φ，使得r（φ，T）＝Pφ－φ＝0.

此外，時(shí)滯微分方程的周期解也由周期T決定.注意到周期解的相位漂移也會產(chǎn)生其他周期解，因此為了獲得精確的周期，需要補(bǔ)充相條件s（φ，T）＝0消除由于相位漂移而產(chǎn)生的不確定性.根據(jù)上述分析，為了求時(shí)滯方程（1）的周期解，需要考慮如下方程：

本文將問題（3）轉(zhuǎn)化為求解如下約束最優(yōu)化問題：尋找初始函數(shù)φ和周期T，在約束條件s（φ，T）＝0下，求函數(shù)‖xT（φ）－φ‖2的最小值，即

為了求解問題（4），定義懲罰函數(shù)：

這里：‖·‖為2－范數(shù)；參數(shù)σ稱為懲罰引子，是一個(gè)較大的正常數(shù).于是，把問題（4）轉(zhuǎn)化為無約束最優(yōu)化問題：

為了尋找快速收斂的無約束最優(yōu)化方法，假設(shè)P＝（φ，T）T是J（φ，T）極小值點(diǎn)的近似，則根據(jù)Newton法，下一次迭代的近似解為

這里：▽2J（P）表示在近似解處的Hessen矩陣；▽J（P）為梯度；而

稱為Newton方向，是下一步迭代的搜索方向.

方程（1）的相空間是無限維的，在時(shí)滯較大的情形下，映射（3）可能會失去緊性.因此，必須在有限維空間［14］下考慮該系統(tǒng).數(shù)值計(jì)算時(shí)滯方程（1）周期解的第一步是對初始函數(shù)φ進(jìn)行離散化.選擇一個(gè)網(wǎng)格步長h＝τ／（N－1），定義網(wǎng)格點(diǎn)τi＝－τ＋（i－1）h.用N 個(gè)網(wǎng)格點(diǎn)上的值φi（i＝1，2，…，N）近似地逼近初始函數(shù)φ.注意離散初始函數(shù)時(shí)也可使用非均勻網(wǎng)格點(diǎn)，不影響計(jì)算效率［15］.

計(jì)算時(shí)滯方程的周期解就是找到合適的初始函數(shù)φ及周期T，在一般打靶法中，通常用函數(shù)插值的方法用N個(gè)φi值逼近初始函數(shù)φ.但當(dāng)周期解不是強(qiáng)烈吸引時(shí)，或者周期解不穩(wěn)定時(shí)，就需要非常多的φi（即N非常大）確保逼近的有效性.而當(dāng)N非常大時(shí)，用Newton法進(jìn)行優(yōu)化計(jì)算時(shí)，需要計(jì)算N×N個(gè)元素的Jacobi矩陣，計(jì)算量龐大.

事實(shí)上，對于給定的初始函數(shù)值φi，需要從這些數(shù)據(jù)中抽取蘊(yùn)含在其中的規(guī)律作為真實(shí)初始函數(shù)的近似，因此通過φi構(gòu)造的初始函數(shù)φ只需能反映φi的變化趨勢即可，不必要求初始函數(shù)在每個(gè)逼近步都精確地通過φi.如圖1所示，本文提出用最小二乘函數(shù)擬合方法構(gòu)造初始函數(shù)φ，這樣構(gòu)造的初始函數(shù)是在能量范數(shù)意義下φi的最佳逼近.

圖1 初始函數(shù)的構(gòu)造方式（虛線為經(jīng)典的（線性）插值方法；實(shí)線為本文提出的擬合方法）Fig.1 Construction of the initial function（the dotted line is the figure of classic（linear）interpolation method and solid lines is the figure of fitting method）

于是，

E的最小值滿足：

進(jìn)而有

本文提出的應(yīng)用最優(yōu)化方法計(jì)算時(shí)滯方程周期解的步驟如下：

選取初始數(shù)據(jù)；給定初始懲罰因子σ＞0，放大系數(shù)β＞1，允許誤差ε1＞0，ε2＞0.

1）對初始函數(shù)φ進(jìn)行離散化.選擇一個(gè)網(wǎng)格步長h＝τ／（N－1），定義τi＝－τ＋（i－1）h.

7）保存初始值和周期，迭代結(jié)束.

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

例1 小世界網(wǎng)絡(luò)的群體動力學(xué)，其運(yùn)動方程［16］如下：

選擇ξ＝3，τ＝1，可以計(jì)算出μ＝0.040 8，T＝4.應(yīng)用本文的最優(yōu)化計(jì)算方法，在區(qū)間［－1，0］隨機(jī)選擇11個(gè)初始函數(shù)節(jié)點(diǎn)，用3階多項(xiàng)式作為擬合基函數(shù)擬合初始函數(shù).系統(tǒng)的周期選為4，因?yàn)橹芷谝阎?，因此相條件s（φ，T）＝0自動滿足.系統(tǒng)的周期解計(jì)算結(jié)果如圖2所示.由圖2可見，在［T－τ，T］上的解片段xT（φ）與初始函數(shù)φ非常接近，相對誤差為

約為0.01%，因此找到了近似周期解.

圖2 用本文方法計(jì)算得到的系統(tǒng)周期解（A）及［T－τ，T］上解xT（φ）和初始函數(shù)的比較（B）Fig.2 Periodic solutions（A）obtained by solving optimization problems and the solution segment at［T－τ，T］and the initial function（B）

對于不同的未知數(shù)個(gè)數(shù)（N）及不同的擬合基函數(shù)（m），表1列出了計(jì)算結(jié)果的相對誤差.由表1可見：對于相同的擬合基函數(shù)m及不同未知數(shù)個(gè)數(shù)N，計(jì)算結(jié)果的誤差幾乎相同；而對于相同的未知數(shù)個(gè)數(shù)N，增加擬合基函數(shù)的項(xiàng)數(shù)m，可以使計(jì)算誤差減小.但對于3階或4階的擬合基函數(shù)，計(jì)算結(jié)果的相對近似約為0.01%，計(jì)算精度較高.因此，可以用低階擬合基函數(shù)及較少的未知數(shù)獲得較好的計(jì)算精度.

表1 不同未知數(shù)個(gè)數(shù)與擬合階的計(jì)算誤差比較（%）Table 1 Relative error for different unknown variables and fitting orders（%）

例2 考慮系統(tǒng)：

根據(jù)文獻(xiàn)［15］，當(dāng)α＝π／2時(shí)，系統(tǒng)（15）產(chǎn)生 Hopf分岔；當(dāng)α＞π／2（并接近π／2）時(shí)，系統(tǒng)（15）具有周期為4的周期解.選擇a＝1.571，應(yīng)用本文的最優(yōu)化計(jì)算方法，在區(qū)間［－1，0］內(nèi)隨機(jī)選擇11個(gè)初始節(jié)點(diǎn)，用4階多項(xiàng)式作為擬合基函數(shù)擬合初始函數(shù)，周期選為4，計(jì)算所得周期解如圖3所示.由圖3可見，在［T－τ，T］上的解片段xT（φ）與初始函數(shù)非常接近，計(jì)算所需的CPU時(shí)間為7.405s，計(jì)算的相對誤差為0.070 2%，因此找到了近似周期解（圖3（A））.如果用插值方法逼近初始函數(shù)，將［－1，0］區(qū)間劃分成11個(gè)節(jié)點(diǎn)（與本文方法未知數(shù)個(gè)數(shù)相同），則系統(tǒng)的計(jì)算誤差僅為0.267 2%，有較大的計(jì)算誤差（圖3（B））.如果要達(dá)到與本文方法相當(dāng)?shù)恼`差（約0.07%），通過計(jì)算可知，至少需將［－1，0］區(qū)間劃分成21個(gè)節(jié)點(diǎn)，計(jì)算所需的CPU時(shí)間為12.788s，計(jì)算量約提升42%.

圖3 本文方法（A）與傳統(tǒng)方法（B）的計(jì)算結(jié)果比較Fig.3 Comparison of computational results calculated by the proposed method（A）and conventional method（B）

例3 考慮具有時(shí)滯的非線性Duffing振子：

選取k＝0.64，α＝0.36，β＝－0.01，計(jì)算可知：當(dāng)τ＞π（并在π附近）時(shí)，系統(tǒng)（16）產(chǎn)生周期為2π的周期解；當(dāng)β＜0，系統(tǒng)（16）產(chǎn)生的周期解是不穩(wěn)定的.表2列出了不同方法的計(jì)算誤差與所用的CPU時(shí)間比較.由表2可見，用本文提出的方法，當(dāng)用3階擬合多項(xiàng)式時(shí)，并不能計(jì)算得出周期解；而當(dāng)m＝5或m＝6時(shí)，可以找到周期解，并且計(jì)算的誤差小于0.01%，表明本文方法對于不穩(wěn)定周期解的計(jì)算有效.如果用傳統(tǒng)的線性插值方法計(jì)算近似周期解，表2的計(jì)算結(jié)果表明，當(dāng)未知數(shù)個(gè)數(shù)為71時(shí)，才能找到合適的初始函數(shù)，花費(fèi)的CPU時(shí)間為157.213s，遠(yuǎn)大于本文方法（N＝15，m＝5）需要的CPU時(shí)間42.658s，計(jì)算時(shí)間提升近72.8%.因此本文方法節(jié)省了計(jì)算量.

表2 不同方法的計(jì)算誤差比較Table 2 Comparison of computational errors of different methods

綜上，本文提出了應(yīng)用最優(yōu)化的方法求解時(shí)滯微分方程的周期解，考慮了周期解的相位漂移條件，將計(jì)算周期解的問題轉(zhuǎn)化為一個(gè)具有約束的最優(yōu)化問題.為了保持Poincaré映射的緊性，求初始函數(shù)的有限維逼近，最小二乘擬合的方法通過有限個(gè)初始函數(shù)值φi逼近真實(shí)的初始函數(shù).數(shù)值實(shí)驗(yàn)表明，用本文的函數(shù)擬合方法，用少量的未知數(shù)，采用適當(dāng)階數(shù)的擬合基函數(shù)即可獲得周期解的初始函數(shù).當(dāng)周期解的吸引性較強(qiáng)時(shí)，采用低階的擬合基函數(shù)即可；當(dāng)周期解的吸引性較弱或周期解不穩(wěn)定時(shí)，可以采用高階的擬合基函數(shù)得到好的逼近效果.由于在本文提出的方法中，未知數(shù)比傳統(tǒng)線性插值方法的未知數(shù)個(gè)數(shù)少很多，因此提高了計(jì)算效率.

［1］XU Xu，WANG Zai－h(huán)ua.Effect of Small World Connection on the Dynamics of a Delayed Ring Network ［J］.Nonlinear Dynamics，2009，56（1／2）：127－144.

［2］WANG Zai－h(huán)ua，HU Hai－yan.Stability and Bifurcation of Delayed Dynamics Systems：From Theory to Applications［J］.Advances in Mechanics，2013，43（1）：3－20.（王在華，胡海巖.時(shí)滯動力系統(tǒng)的穩(wěn)定性與分岔：從理論走向應(yīng)用 ［J］.力學(xué)進(jìn)展，2013，43（1）：3－20.）

［3］Kuang K.Delay Differential Equations：With Applications in Population Dynamics［M］.New York：Academic Press，1993.

［4］Luzyanina T，Engelborghs K，Lust K，et al.Computation，Continuation and Bifurcation Analysis of Periodic Solutions of Delay Differential Equations ［J］.International Journal of Bifurcation and Chaos，1997，7（11）：2547－2560.

［5］Lust K，Spence A，Champneys A R.An Adaptive Newton－Picard Algorithm with Subspace Iteration for Computing Periodic Solutions［J］.SIAM Journal Scientific Computations，1998，19（4）：1188－1209.

［6］HE Ying－sheng.Numerical Solution and Matlab Simulation of Impulsive Delay Differential Equations［J］.Journal of Jishou University：Natural Sciences Edition，2009，30（4）：30－33.（何迎生.脈沖時(shí)滯微分方程的數(shù)值解法及Matlab實(shí)現(xiàn) ［J］.吉首大學(xué)學(xué)報(bào)：自然科學(xué)學(xué)版，2009，30（4）：30－33.）

［7］Minamoto T，Nakao M T.A Numerical Verification Method for a Periodic Solution of a Delay Differential Equation［J］.Journal of Computational and Applied Mathematics，2010，235（3）：870－878.

［8］WANG Wan－yong，XU Jian.Multiple Scales Analysis for Double Hopf Bifurcation with 1∶3Resonances［J］.Nonlinear Dynamics，2011，66（1／2）：39－51.

［9］YOU Xiang－chang，XU Hang.Analytical Approximations for the Periodic Motion of the Duffing System with Delayed Feedback［J］.Numerical Algorithms，2011，56（4）：561－576.

［10］WANG Wan－yong.Resonant Double Hopf Bifurcation and Its Singularity in Systems with Delay Coupling ［D］.Shanghai：Tongji University，2011.（王萬永.時(shí)滯耦合系統(tǒng)共振雙 Hopf分岔奇異性及其分類 ［D］.上海：同濟(jì)大學(xué)，2011.）

［11］ZHENG Yuan－guang.Nonlinear Dynamics of Singularly Perturbed Dynamical Systems with Delays［D］.Nanjing：Nanjing University of Aeronautics and Astronautics，2009.（鄭遠(yuǎn)廣.含時(shí)滯的奇異攝動系統(tǒng)的非線性動力學(xué)［D］.南京：南京航空航天大學(xué)，2009.）

［12］SHENG Xia，ZHANG Li，ZHENG Hong，et al.Periodic Solutions of a Higher Dimensional Network with Time Delay［J］.Journal of Jilin University：Science Edition，2010，48（5）：728－732.（盛夏，張麗，鄭虹，等.具有時(shí)滯的高維網(wǎng)絡(luò)系統(tǒng)的周期解 ［J］.吉林大學(xué)學(xué)報(bào)：理學(xué)版，2010，48（5）：728－732.）

［13］FENG Zi－xuan，JI Shu－guan，XU Xu.Finding Periodic Solution of Differential Equation via Solving Optimization［J］.Journal of Jilin University：Science Edition，2009，47（1）：37－39.（馮子玹，冀書關(guān)，徐旭.通過求解最優(yōu)化問題尋找微分方程的周期解 ［J］.吉林大學(xué)學(xué)報(bào)：理學(xué)版，2009，47（1）：37－39.）

［14］Küpper T，LI Yong，ZHANG Bo.Periodic Solutions for Dissipative－Repulsive Systems ［J］.Tohoku Mathematical Journal，2000，52（3）：321－329.

［15］XU Jian，Chung K W.Dynamics for a Class of Nonlinear Systems with Time Delay［J］.Chaos，Solitons ＆Fractals，2009，40（1）：28－49.

［16］LI Chun－gang，CHEN Guan－rong.Local Stability and Hopf Bifurcation in Small－World Delayed Networks［J］.Chaos，Solutons ＆ Fractals，2004，20（2）：353－361.