賈 蒙

(新鄉(xiāng)學(xué)院 機(jī)電工程學(xué)院,河南 新鄉(xiāng) 453003)

高維非線性映射系統(tǒng)的不穩(wěn)定流形計(jì)算方法研究

賈 蒙

(新鄉(xiāng)學(xué)院 機(jī)電工程學(xué)院,河南 新鄉(xiāng) 453003)

提出了一種基于曲率約束條件的計(jì)算離散動(dòng)力系統(tǒng)鞍型不動(dòng)點(diǎn)一維不穩(wěn)定流形的新算法，并以Hénon映射為例進(jìn)行了計(jì)算。新算法以增長(zhǎng)流形為基本思想，通過曲率約束和距離控制來(lái)確定離散點(diǎn)間的距離；提出流形的偏轉(zhuǎn)角度可以通過流形上的已知點(diǎn)來(lái)預(yù)測(cè)，解決了流形上新點(diǎn)的原像位置快速確定的困難。仿真發(fā)現(xiàn)：Hénon映射的一維不穩(wěn)定流形在標(biāo)準(zhǔn)參數(shù)下與Hénon映射產(chǎn)生的散點(diǎn)圖分布一致，在其它幾組參數(shù)下，一維不穩(wěn)定流形的兩個(gè)分支之間保持著某種程度的對(duì)稱性，該研究對(duì)Hénon映射的進(jìn)一步研究打下基礎(chǔ)。

離散動(dòng)力系統(tǒng);雙曲不動(dòng)點(diǎn);不穩(wěn)定流形;Hénon映射;混沌

動(dòng)力系統(tǒng)按照其狀態(tài)變量隨時(shí)間的變化是否連續(xù)可分為兩類：連續(xù)動(dòng)力系統(tǒng)和離散映射動(dòng)力系統(tǒng)。其中離散動(dòng)力系統(tǒng)可看作是對(duì)連續(xù)動(dòng)力系統(tǒng)采樣后的結(jié)果，而離散動(dòng)力系統(tǒng)經(jīng)常以迭代函數(shù)的形式出現(xiàn)?？茖W(xué)研究及生活實(shí)際中的很多過程都可以用迭代函數(shù)的形式來(lái)描述，Hénon映射[1]就是一個(gè)比較著名的例子，它具有一定的混沌行為，在圖像加密[2]、分形研究[3]等領(lǐng)域都有廣泛地應(yīng)用。其表達(dá)式為

(1)

將式(1)寫成映射函數(shù)的形式

(2)

無(wú)論是連續(xù)時(shí)間系統(tǒng)，還是離散時(shí)間系統(tǒng)，雙曲平衡點(diǎn)或雙曲周期軌道的穩(wěn)定、不穩(wěn)定流形的計(jì)算對(duì)動(dòng)力學(xué)問題的全局性理解起到了非常重要的作用。 對(duì)于一個(gè)給定的非線性動(dòng)力學(xué)問題，進(jìn)行嚴(yán)格的理論分析是非常復(fù)雜和艱難的，通常采用數(shù)值方法計(jì)算其穩(wěn)定和不穩(wěn)定流形。

Krauskopf等[4]提出利用測(cè)地線距離作為參數(shù)，不需要對(duì)向量場(chǎng)進(jìn)行調(diào)整，而是通過求解一系列邊值問題以計(jì)算流形。 它能夠很好的控制流形計(jì)算中的增長(zhǎng)問題，但是邊值問題的求解計(jì)算量很大，特別當(dāng)與向量場(chǎng)近似相切的時(shí)候所需要的步驟會(huì)更多，耗時(shí)嚴(yán)重． 另外，該方法可以以徑向的軌線構(gòu)建流形，但不能繪制出動(dòng)力學(xué)系統(tǒng)的流． 其他方法，如 PDE 方法[5]、軌道參數(shù)化方法[6]、細(xì)分法[7]等，在計(jì)算流形過程中沒有給出有效控制流形各個(gè)方向增長(zhǎng)的方法。

本文提出了一種基于曲率約束條件的計(jì)算離散動(dòng)力系統(tǒng)鞍型不動(dòng)點(diǎn)一維不穩(wěn)定流形的新算法，解決了流形上新點(diǎn)的原像位置快速確定的困難。

1 流形計(jì)算的基本方法

首先我們將研究范圍擴(kuò)展到一般映射。設(shè)F:Rn→Rn為一個(gè)保向的微分同胚映射函數(shù)，

xm+1=F(xm)

(3)

取x0為映射的初始點(diǎn)，則系統(tǒng)滿足

Fk+l(x0)=Fk(Fl(x0))

(4)

式中:k,l為任意正整數(shù);當(dāng)F可逆時(shí);k,l可以取任意整數(shù)。

對(duì)于映射式(3)，若存在x0∈Rn，滿足F(x0)=x0，則稱x0為系統(tǒng)的不動(dòng)點(diǎn)。在x0處對(duì)系統(tǒng)進(jìn)行線性化，將式(3)轉(zhuǎn)化為線性映射x→Ax，其中A為F在x0處的雅各比矩陣A=DF(x0)=[?fi/?xj](x0)。若矩陣A的特征值的模都不等于1，那么x0就是一個(gè)雙曲不動(dòng)點(diǎn)。其中模小于1的特征值叫做穩(wěn)定特征值，其對(duì)應(yīng)的特征向量{v1,v2,…,vl}張成穩(wěn)定特征空間Es；模大于1的特征值叫做不穩(wěn)定特征值，它們對(duì)應(yīng)的特征向量{vl+1,vl+2,…,vn}張成不穩(wěn)定特征空間Eu。全空間

Rn=Es⊕Eu

(5)

設(shè)x0是微分同胚映射函數(shù)F的一個(gè)雙曲不動(dòng)點(diǎn)，則在x0的鄰域U內(nèi)存在局部穩(wěn)定和不穩(wěn)定流形：

(6)

(7)

定理的證明見文獻(xiàn)[4]。全局穩(wěn)定和不穩(wěn)定流形定義為

穩(wěn)定流形和不穩(wěn)定流形在分析系統(tǒng)的動(dòng)力學(xué)特性中起著非常重要的作用，它們充當(dāng)不同吸引子吸引域的邊界，將全空間劃分為多個(gè)具有不同動(dòng)力特性的不變子空間，而當(dāng)穩(wěn)定與不穩(wěn)定流形相交時(shí)，就會(huì)引起同宿、異宿以及混沌[3-6]等復(fù)雜動(dòng)力學(xué)行為的出現(xiàn)。不變流形計(jì)算對(duì)于上述問題的研究都有著積極地促進(jìn)作用。

離散映射動(dòng)力系統(tǒng)的流形計(jì)算比較困難，對(duì)于一維流形的計(jì)算，比較著名的算法見文獻(xiàn)[7-8]。Krauskopf等[9-10]通過每步向流形中加入一個(gè)離散點(diǎn)來(lái)增長(zhǎng)流形，離散點(diǎn)間的距離受流形局部曲率的控制，該算法能夠有效展示流形的細(xì)節(jié)。該算法的不足之處在于，它用二分法搜索已有流形區(qū)間段來(lái)確定新離散點(diǎn)的前像位置，比較繁瑣，這也是影響該算法計(jì)算速度的一個(gè)關(guān)鍵因素。

本文提出的算法從增長(zhǎng)流形[11-12]這個(gè)思想出發(fā)的，提出了一種預(yù)測(cè)流形偏轉(zhuǎn)角度的方法，新方法可以快速確定流形上新離散點(diǎn)的原像位置，同時(shí)采用了中曲率控制[13-14]和距離控制條件，算法簡(jiǎn)潔明了，計(jì)算速度快，且易于編程實(shí)現(xiàn)。

2 一維流形計(jì)算方法的誤差分析

賈蒙介紹的各種一維流形計(jì)算方法，都沒有對(duì)算法的計(jì)算精度進(jìn)行討論。因?yàn)橛?jì)算映射的一維流形時(shí)，經(jīng)常會(huì)遇到混沌映射，而混沌的一個(gè)顯著特點(diǎn)就是對(duì)于初始誤差的敏感性，意味著計(jì)算誤差將隨計(jì)算的進(jìn)行而指數(shù)放大，對(duì)于研究流形計(jì)算方法的人來(lái)說，這不啻為一場(chǎng)災(zāi)難。

既然誤差會(huì)指數(shù)放大，所以采用精度有限、具有截?cái)嗾`差的計(jì)算機(jī)來(lái)，計(jì)算流形似乎是不可能的。但其實(shí)沒有這么絕望?；煦缦到y(tǒng)雖然對(duì)于初值具有敏感依賴性，長(zhǎng)期行為近似隨機(jī)，但短期行為還是可以預(yù)測(cè)的，也就是說，對(duì)長(zhǎng)期行為進(jìn)行計(jì)算的話誤差較大，但計(jì)算短期行為時(shí)誤差仍然是保持有界的。而我們?cè)谟?jì)算映射系統(tǒng)的一維流形時(shí)，計(jì)算的弧長(zhǎng)都非常有限，因此迭代的次數(shù)也很少，所以其誤差能夠保持在有界范圍內(nèi)。

并且，一切形式為

L={(x,y):|x|<δ,y=p(x)}

的流形L的迭代都將在C1拓?fù)湎率諗坑诓环€(wěn)定流形，即

(10)

設(shè)M是Rn空間中的一個(gè)m維子流形，在M的每個(gè)點(diǎn)x∈M處定義維數(shù)為n-m的超平面N(x)，該超平面在x處與M橫截相交。不失一般性，假設(shè)N(x)是光滑的，而且N(x)的選擇并不唯一。集合E={N(x),x∈M}在M上具有向量纖維叢(vector fiber bundle)結(jié)構(gòu)，因此可以把M看作零截面(null section){(x,0)}。定義(x⊕y)=x+y，其中x∈M、y∈N(x)。對(duì)于任意固定的小量ε>0，稱集合Uε={(x⊕y):x∈M,y∈N(x),|y|<ε}為子流形M的一個(gè)管狀鄰域。

設(shè)Nε是流形M的一個(gè)管狀鄰域，流形M1?Nε的形式為

M1={x⊕h(x):h(x)∈N(x)}

其中x⊕h(x):M→Rn是一個(gè)光滑映射。滿足該條件的M1被稱作是處于M的C1鄰域內(nèi)。C1距離定義為映射的C1范數(shù)(C1norm)

(11)

令原點(diǎn)0是微分同胚F:Rn→Rn的一個(gè)雙曲不動(dòng)點(diǎn)，則任意形式為

L={(x,y):|x|<δ,y=p(x)}

的流形L將C1拓?fù)涫諗坑诓环€(wěn)定流形Wu的緊致部分，即當(dāng)m→∞時(shí)，Wm=Fk(L)∩Nε→M。迭代將逐點(diǎn)收斂于全局不穩(wěn)定流形Wu。

當(dāng)用逆函數(shù)F-1替換函數(shù)F時(shí)，以上理論也可適用于穩(wěn)定流形。

3 基于曲率約束條件的流形計(jì)算方法

考慮軌道上相鄰三個(gè)離散點(diǎn)之間的角度，如圖2所示，由圖易得

(12)

其中

(13)

(14)

檢查α是否滿足下列約束條件

αmin<α<αmax

(Δα)min<Δkα<(Δα)max

(15)

將pk+1加入到序列M中，并判斷M′的長(zhǎng)度arcl是否已經(jīng)達(dá)到所需弧長(zhǎng)ARC，滿足則結(jié)束，不滿足則繼續(xù)計(jì)算。

上述步驟都結(jié)束后，再沿A的不穩(wěn)定特征向量的反方向計(jì)算Wu(x0)的另一分支，兩個(gè)分支之和就是Wu(x0)。

4 高維不穩(wěn)定流形計(jì)算與仿真

4.1三維Hénon映射的一維不穩(wěn)定流形

三維Hénon映射的表達(dá)式為

(16)

當(dāng)M1=1.4,M2=0.2,B=0.1時(shí)，具有與二維Hénon映射類似的混沌吸引子，如圖2(a)所示。易得該映射具有不動(dòng)點(diǎn)

x0=(x*,y*,z*)，其中

該不動(dòng)點(diǎn)是二次映射F2的一個(gè)雙曲不動(dòng)點(diǎn)，且具有一維不穩(wěn)定流形，計(jì)算結(jié)果如圖2(b)所示。

(a) 混沌吸引子

(b) 一維不穩(wěn)定流形，弧長(zhǎng)為30，計(jì)算耗時(shí)0.6 s

觀察比較圖2中的兩幅圖，發(fā)現(xiàn)二者幾乎一模一

樣，不同之處在于：(a)圖中的點(diǎn)是混沌的，近似隨機(jī)排列的，如果按照迭代的順序?qū)⑦@些點(diǎn)順序連接起來(lái)的話，整個(gè)圖將是雜亂無(wú)章的，而(b)圖中的點(diǎn)是有序的。我們發(fā)現(xiàn)這點(diǎn)同樣適用于二維Hénon映射。

當(dāng)M1=0.19、M2=0.999 1、B=0時(shí),x0也是一個(gè)雙曲不動(dòng)點(diǎn)，一維流形計(jì)算結(jié)果如圖3所示。

圖3 三維Hénon映射的一維不穩(wěn)定流形

Fig.3 One dimensional unstable manifold of three dimensional Hénon system

圖3中間似乎“雜亂”卻又左右對(duì)稱的部分為三維Hénon映射的吸引子，可見該吸引子完全被一維不穩(wěn)定流形所“包圍”。

4.2四維映射的一維流形計(jì)算

為了說明新算法能夠用于高維映射的流形計(jì)算，特將三維Hénon映射增加了一維，變?yōu)樗木SHénon映射，表達(dá)式為

(17)

(a) w=0

(b) z=0

(c) y=0

(d) x=0

圖4 四維Hénon map的一維不穩(wěn)定流形

Fig.4 One dimensional unstable manifold of four dimensional Hénon system

5 總 結(jié)

本文提出了一種計(jì)算離散動(dòng)力系統(tǒng)不動(dòng)點(diǎn)一維不穩(wěn)定流形的新算法，該算法通過每步加入一個(gè)離散點(diǎn)來(lái)增長(zhǎng)流形，相鄰離散點(diǎn)間的距離通過曲率約束和距離控制條件來(lái)確定。提出了一種不穩(wěn)定流形上角度的預(yù)測(cè)方法，可以用來(lái)快速地確定當(dāng)前所要加入點(diǎn)的原像的位置，這是相對(duì)于文獻(xiàn)[9]中提出的算法的優(yōu)越之處。此外，本算法可同時(shí)用于一維穩(wěn)定和不穩(wěn)定流形的計(jì)算，并且導(dǎo)數(shù)傳遞的推導(dǎo)在高維空間下依舊成立，算例5.1和5.2說明了本章算法可以用于高維空間下的一維流形計(jì)算。在計(jì)算一維穩(wěn)定流形時(shí)，本章算法不需要映射的逆函數(shù)F-1的顯式表達(dá)式，所以通用性也較傳統(tǒng)算法有了提高。

本算法需要計(jì)算映射函數(shù)的Jacobian矩陣，所以得考慮該矩陣的計(jì)算難度。對(duì)于顯式定義的同胚映射函數(shù)，其Jacobian矩陣可以寫成顯式表達(dá)式，對(duì)于某個(gè)特定點(diǎn)，只要將坐標(biāo)帶入即可求得，計(jì)算比較簡(jiǎn)單，耗時(shí)與函數(shù)進(jìn)行一次迭代相當(dāng)；但對(duì)于向量場(chǎng)的Poincaré映射，其Jacobian矩陣需要用數(shù)值方法進(jìn)行求解，比較復(fù)雜，在這種情況下，本算法可能是不具優(yōu)勢(shì)的。對(duì)于非同胚映射，或者函數(shù)不可逆以及有多個(gè)逆，本算法都不能夠使用。

另外，當(dāng)公式(1)參數(shù)a和b取其它更多變化時(shí)，Hénon映射的不穩(wěn)定流形可能會(huì)表現(xiàn)出更加豐富的性質(zhì)以及更加漂亮的圖形，這有待我們進(jìn)一步的探索驗(yàn)證。我們希望本文的計(jì)算結(jié)果能對(duì)Hénon映射的進(jìn)一步研究產(chǎn)生一些啟發(fā)。

[1] 賈蒙. 映射動(dòng)力系統(tǒng)一維流形并行計(jì)算方法[J]. 振動(dòng)與沖擊, 2014,33(9):40-47.

JIA Meng. A parallel algorithm for approximating 1-D manifold of map[J].Journal of Vibration and Shock, 2014,33(9):40-47.

[2] GAO H, ZHANG Y, LIANG S, et al. New chaotic algorithm for image encryption[J].Chaos Solutions and Fractals,2006(29): 393-399.

[3] GRASSBERGER P. On the fractal dimension of the Henon attractor[J]. Physics Letters A, 1983, 97(6):224-226.

[4] PALIS J, MELO W D. Geometric theory of dynamical Systems[M].Springer, 1982.

[5] YOU Z, KOSTELICH E J, YORKE J A. Calculating stable and unstable manifolds[J]. Int. J. Bifurc. Chaos Appl.Sci. Eng.,1991, 1(3): 605-623.

[7] PARKER T S, CHUA L O. Practical numerical algorithms for Chaotic Systems[M]. Springer, Berlin, 1989.

[8] HOBSON D. An efficient method for computing invariant manifolds[J]. J. Comput. Phys. 1991,104: 14-22.

[9] KRAUSKOPF B, OSINGA H M. Growing unstable manifolds of planar maps[J]. 1517, 1997, http://www.ima.umn.edu/preprints/OCT97/1517.ps.gz.

[10] KRAUSKOPF B, OSINGA H M. Growing 1D and quasi-2D unstable manifolds of maps[J]. J. Comput. Phys.1998,146: 406-419.

[11] ENGLAND J P, KRAUSKOPF B, OSINGA H M. Computing One-Dimensional Stable Manifolds and Stable Sets of Planar Maps without the Inverse[J]. SIAM J. Appl. Dyn. Syst. 2004,3(2):161-190.

[12] DELLNITZ M, HOHMANN A. A subdivision algorithm for the computation of unstable manifolds and global attractors[J]. Numer. Math, 1997,75: 293-317.

[13] FUNDINGER D. Toward the calculation of higher-dimensional stable manifolds and stable sets for noninvertible and piecewise-smooth maps[J]. J Nonlinear Sci, 2008, 18: 391-413.

[14] KRAUSKOPF B, OSINGA H M. Globalizing two-dimensional unstable manifolds of maps[J]. Int. J. Bifurc.Chaos,1998, 8(3): 483-503.

Computationmethodforunstablemanifoldofhénonmap

JIA Meng

(College of Mechanical & Electrical Engineering， Xinxiang College， Xinxiang 453003, China)

A new algorithm was presented for computing one dimensional unstable manifold for saddle fixed point of a discrete dynamic system based on the constraint condition of curvature. Hénon map was taken as an example to check the performance of the algorithm. The manifold growing was taken as the basic idea of the new algorithm. Curvature constraint and distance control were used to determine the distance between discrete points. The unstable manifold grew with new point added at each step and the distance between consecutive points was adjusted according to the local curvature. It was proved that the gradient of the manifold at the new point can be predicted with the known points on the manifold and in this way the preimage of the new point can be located immediately. The simulation showed that the one dimensional unstable manifold of Hénon map coincides with the scatter point diagram distribution produced by itself under standard parameters; under the other several groups of parameters, two branches of the unstable manifold are nearly symmetric, and they serve as the borderline of Hénon map iteration sequence. The results laid a foundation for further studying Hénon map.

discrete dynamic system; hyperbolic fixed point; unstable manifold; Hénon map; chaos

國(guó)家自然科學(xué)基金(61501391);河南省高等學(xué)校重點(diǎn)科技項(xiàng)目(15A510035)

2016-03-14 修改稿收到日期：2016-07-05

TP301.6

： A

10.13465/j.cnki.jvs.2017.17.038

作 者 賈蒙 男，博士，副教授，1981年12月生