賈 蒙

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

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

賈 蒙

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

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

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

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

(1)

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

(2)

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

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

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

1 流形計算的基本方法

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

xm+1=F(xm)

(3)

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

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

(4)

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

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

Rn=Es⊕Eu

(5)

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

(6)

(7)

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

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

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

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

2 一維流形計算方法的誤差分析

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

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

并且，一切形式為

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

的流形L的迭代都將在C1拓撲下收斂于不穩(wěn)定流形，即

(10)

設(shè)M是Rn空間中的一個m維子流形，在M的每個點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)。對于任意固定的小量ε>0，稱集合Uε={(x⊕y):x∈M,y∈N(x),|y|<ε}為子流形M的一個管狀鄰域。

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

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

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

(11)

令原點0是微分同胚F:Rn→Rn的一個雙曲不動點，則任意形式為

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

的流形L將C1拓撲收斂于不穩(wěn)定流形Wu的緊致部分，即當m→∞時，Wm=Fk(L)∩Nε→M。迭代將逐點收斂于全局不穩(wěn)定流形Wu。

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

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

考慮軌道上相鄰三個離散點之間的角度，如圖2所示，由圖易得

(12)

其中

(13)

(14)

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

αmin<α<αmax

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

(15)

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

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

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

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

三維Hénon映射的表達式為

(16)

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

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

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

(a) 混沌吸引子

(b) 一維不穩(wěn)定流形，弧長為30，計算耗時0.6 s

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

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

當M1=0.19、M2=0.999 1、B=0時,x0也是一個雙曲不動點，一維流形計算結(jié)果如圖3所示。

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

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

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

4.2四維映射的一維流形計算

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

(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é)

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

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

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

[1] 賈蒙. 映射動力系統(tǒng)一維流形并行計算方法[J]. 振動與沖擊, 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

國家自然科學(xué)基金(61501391);河南省高等學(xué)校重點科技項目(15A510035)

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

TP301.6

： A

10.13465/j.cnki.jvs.2017.17.038

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