基于Foliation條件的離散動力系統二維流形計算

賈 蒙

（新鄉學院機電工程學院，河南新鄉 453003）

基于Foliation條件的離散動力系統二維流形計算

賈 蒙

（新鄉學院機電工程學院，河南新鄉 453003）

研究離散動力系統雙曲不動點的二維流形計算，利用不變流形軌道上Jacobian矩陣能夠傳遞導數這一特殊性質，提出一種新的一維流形計算方法，通過預測－校正兩個步驟迅速確定流形上新網格點，避免重復計算，并簡化精度控制條件.在此基礎上，將基于流形面Foliation條件進行推廣，推廣后的Foliation條件能夠控制二維流形上的一維子流形的增長速度，從而實現二維流形在各個方向上的均勻增長.此外，算法可以同時用于二維穩定和不穩定流形的計算.以超混沌三維Hénon映射和具有蝶形吸引子的Lorenz系統為例驗證了算法的有效性.

離散動力系統；穩定流形；不穩定流形；導數傳遞；三維Hénon映射；Lorenz系統；混沌吸引子

0 引言

穩定流形和不穩定流形在分析系統的動力學特性中起著非常重要的作用，它們充當著不同吸引子吸引域的邊界，將全空間劃分為多個具有不同動力特性的不變子空間；而當穩定與不穩定流形相交時，就會引起同宿、異宿以及混沌等復雜動力學行為的出現.

流形計算作為動力系統分析的一種幾何方法，一直以來都受到研究人員的重視.研究的對象主要有兩類系統：連續動力系統（向量場）和離散動力系統.連續動力系統一般以微分方程的形式來表示，其一維流形只需通過積分求解初值問題就能得到，比較簡單，二維流形計算就變得困難多了：其中一個最主要的問題是當系統的兩個穩定（或不穩定）特征值的比值較大時，流形軌道在各個方向上的增長速度不均，而且軌道朝其中一個特征值方向匯集，從而使得計算的結果變差.針對這個問題，國內外的學者給出了許多不同解決辦法，詳見文獻［1－7］.

離散動力系統一般以映射的形式出現，所以有時也稱作映射動力系統.離散動力系統的流形計算比較困難，與連續動力系統的流形計算相比，它的一個重要特點是軌道上新網格點對軌道上已有的網格點的依賴性，也就是說，新網格點必須是通過已有網格點或其插值經過映射后才能得到；而連續動力系統的流形計算是通過求解初值問題來計算的，所以只有給定一個出發點，就能計算之后的軌道，而無需考慮出發點之前軌道的狀態.這種依賴性的存在，使得離散動力系統的流形計算更為困難.

You［8］，Simó［9］，Parker和Chua［10］，Hobson［11］給出了基域迭代的改進算法.Bernd Krauskopf和Hinke Osinga［12－13］通過每步向流形中加入一個離散點來增長流形，結合曲率控制技術，使得流形的增長速度由流形的局部曲率來決定，該算法能夠有效展示流形的細節.與此思想類似，文獻［14］中提出的算法能夠在逆映射不能顯式表達的情況下計算系統的一維穩定流形，而之前提到的大部分算法都是把系統的穩定流形當作其逆映射的不穩定流形來計算的.Dellnitz和Hohmann提出的細分法［15］也能有效地計算一維流形；Danny Fundinger［16］的算法與此思想類似.以上算法在計算的快速性上都有待提高.一維流形計算的問題沒有有效地解決，所以二維流形的計算就顯得更為困難了，除了要面對計算向量場二維流形時的困難，還要牽扯到插值點原像的搜索問題，并且插值點之后的軌道計算對于插值點之前的軌道存在依賴性，所以就不得不對插值點之前原本不存在的軌道也進行插值，從而使計算變得尤為復雜.目前對于離散系統二維流形的算法很少，文獻［17］中提出了一種應用Foliation條件的算法，該算法比較好地解決了流形各向增長速度不均的問題，實際上與向量場二維流形算法 是同一算法，該算法默認新網格點的原像都處于當前流形邊界上，很多時候這與實際情況是不符的，所以應用受到限制.國內在相關方面的研究較少，文獻［7］提出對于一維流形的計算，最基本的方法是對基域內的有限個離散點進行迭代，這種算法操作簡單，計算速度快，但在迭代過程中點的分離可能導致在基域上相距很近的兩個點經迭代后得到的像間距變得很大，從而錯過流形的細節，使得計算的精度變差.文獻［18］提出一種預估流形增長速度的一維流形算法，在此基礎上，提出了一種快速的二維流形計算方法，但是該算法也沒有解決流形各向增長速度不均的問題進行處理.文獻［21］提出的方法仍然不能反映一維子流形的細節特性.文獻［22］提出了典型Foppl-von Karman（F-K）理論結合全局分形理論進行非線性系統的流形計算，針對某些特定系統取得一定成果，但是該方法仍不能普及到所有的非線性系統流形計算中.

一維流形計算是二維流形計算的基礎，本文在提出一種快速的一維流形算法的基礎上，對Foliation條件進行推廣，提出一種新的二維流形計算方法，該算法很好地保證了二維流形在各個方向上的均勻增長，并巧妙地避免了新網格點的原像搜索問題.

1 流形的基本概念

設F：Rn→Rn為一個保向的微分同胚映射函數

對于系統（1），若存在x0，滿足F（x0） ＝x0，則稱x0為系統的不動點．A為x0處的雅各比矩陣A ＝DF（x0）＝［?fi／?xj］（x0）．若矩陣A的特征值的模都不等于1，那么x0就是一個雙曲不動點；其中模小于1的特征值叫做穩定特征值，其對應的特征向量｛v1，v2，…，vl｝張成穩定特征空間Es；模大于1的特征值叫做不穩定特征值，它們對應的特征向量｛vl＋1，vl＋2，…，vn｝張成不穩定特征空間Eu．全空間Rn＝Es⊕Eu．

定理1設x0是保向微分同胚映射函數F的一個雙曲不動點，則在x0的鄰域U內存在局部穩定和不穩定流形

從上面的定義可知，全局流形就是局部流形在映射F迭代下得到的像．其中全局穩定流形是局部穩定流形在映射F下進行反向迭代得到的像，而全局不穩定流形是局部不穩定流形在映射F下進行正向迭代得到的像，迭代的方向與其直接定義有所區別.

2 一維流形的計算

一維流形簡單看來就是一條滿足特殊條件的曲線，其特殊性之一就在于：對于處于流形軌道上的任意一點，其像和原像必定也處于流形軌道上.本文的算法就是從這一特殊性出發的.

2.1 導數傳遞

由式（4）、（5）易知，對于任意x∈Ws（x0）（或x∈Wu（x0）），有Fk（x）∈Ws（x0）（或Fk（x）∈Wu（x0））其中k可取任意整數．所以我們可以把穩定流形和不穩定流形的計算當作同一類不變集W來研究.

對于任意x∈W（x0），有F（x）∈W（x0）．對式（1）中的函數F在x處進行Taylor級數展開

其中A為F在x處的Jacobian矩陣.對于任意n維系統，有

二維狀態空間下的情形比較容易理解：Δx2／Δx1可看作是一維流形軌道在x處的斜率，經Jacobian矩陣A變換后得到流形軌道在F（x）處的斜率Δx′2／Δx′1，這就是不變集軌道上的導數傳遞，高維狀態空間下的情況類似［23－24］.

在計算穩定流形的過程中，經常用到的是F－1，為了應用的方便，對F在F－1（x）處進行Taylor級數展開，可得

2.2 一維流形計算

以一維不穩定流形的計算為例進行說明.

1）初始化：以不動點x0為起點，沿Eu方向距離為δ處取一點x1，對x1進行映射得到下一點x2＝F（x1），網格點的初始序列為M＝｛x0，x1，x2｝.

2）增長流形.假定當前網格點序列M＝｛x0，…，xn｝，xn＝F（x′n）的原像為x′n，如圖1所示，現要在距xn為Δn增加下一點xn＋1，可以通過以下兩步來完成：

圖1 增長流形Fig.1 Manifold extending

b）校正：Δ′n和Δn之間的比值為R＝Δn／Δ′n，然后對ε進行調整，ε＝Rε．由于ε和ε′是線性關系，所以調整之后ε′的模必然為Δn.

xn＋1＝F（x′n＋ε）．之后要檢查圖1中的夾角α是否滿足以下條件

當條件（13）滿足時，說明步長Δk合適，則接受xn＋1作為流形上的新點．當α＜αmin時，說明步長過小，但仍然接受xn＋1作為流形軌道上的新點，并設定Δn＋1＝SΔn；當α＞αmax時，則說明步長過大，設定Δn＝Δn／S并重新計算xn＋1．其中S為增長因子，為了保證網格點密度不至于變化太快，建議S的取值范圍為1.2＜S＜1.6．此外，當Δn＜Δmin時，為了避免步長過小而影響計算速度，即使條件（13）不滿足，也接受xn＋1．與此相對應，當Δn＞Δmax時，則用線性插值插入一個離散點．當流形弧長達到ARC時計算結束.

一維穩定流形的計算方法與此類似，但是由于牽扯到函數的求逆運算，所以還是存在一些差別.現將二者的不同之處列于下表1.

表1 一維穩定和不穩定流形的計算差異Table 1 Differences between algorithm in com puting 1D unstable and stablemanifolds

3 二維流形的計算

二維流形計算在一定意義上就是一族一維流形軌道計算的集合，有了快速的一維流形算法，二維流形的計算就水到渠成了.現在的主要問題是如何確保二維流形在各個方向上的均勻增長.

3.1 推廣的Foliation條件

Foliation條件在文獻［17］中率先提出，內容如下：在線性Foliation的每個葉形中，與不穩定流形存在唯一的相交線.也就是說，不穩定流形和每個葉形橫截相交.

文獻［17］應用該條件來確保二維流形在各個方向上的均勻增長.每個葉形定義為與流形環面邊界橫截相交的一個超平面，該超平面與流形面存在唯一交線.通過選取交線上與流形邊界特定距離處的網格點來構成下一個流形環面，二維流形以環面的形式向外增長.Foliation條件通過使流形環面每步向外延伸相同的長度來實現二維流形在各個方向上的均勻增長，由于二維流形與一個平面是拓撲等價的，所以如果將二維流形“展平”，則可以看出，Foliation條件是利用畫圓的思想來保證二維流形在各個方向的均勻增長，這個圓的中心為不動點，而圓的半徑每步向外增加相同長度.

Foliation條件通過控制環面的擴張來計算二維流形，而二維流形在本質上是由一維流形軌道“鋪成”的，所以在實際計算過程中就會產生很多的不便.下面我們將從控制流形軌道的增長速度入手，提出推廣的Foliation條件.

首先說明流形軌道上葉形的定義方法：如圖2所示，雙曲不動點x0的二維局部（不）穩定流形面由兩個（不）穩定特征值對應的特征向量V1和V2張成，n＝V1×V2為局部流形面的法向量，對于流形軌道上任意一點xk，由向量和n確定了過點xk的平面，即為該點處的葉形.只要該平面取得足夠小，與流形面的交線必定唯一.

推廣的Foliation條件內容如下：取流形軌道上任意一點，對該點處的切向分量向該點處的葉形進行投影，沿軌道對投影值進行積分，所得結果就是該流形軌道相對于不動點x0的Foliation弧長.若對二維流形面上的所有流形軌道都計算同樣Foliation弧長，就能夠保證二維流形在各個方向上的均勻增長.

由于計算過程中流形軌道通常都是用一系列的離散點來表示，是分段線性的，所以Foliation弧長的積分可以用下面的公式來表示.

圖2 推廣的Foliation條件中葉形的定義Fig.2 Definition of linear foliation in generalized Foliation condition

3.2 二維流形計算步驟

以二維不穩定流形的計算為例.由于二維流形實際上就是無窮多條一維子流形的集合，一維流形的計算是二維流形計算的基礎，不同于單純的一維流形計算，在二維流形中計算一維流形軌道時要按照式（14）標出軌道上離散點的Foliation弧長，從而確保二維流形在各個方向上的均勻增長.以不動點x0為圓心，在局部不穩定流形上以δ為半徑作圓，在圓上均勻取N個點．然后任取其中一點，與x0相連構成初始流形段，然后按照2.2節所介紹的一維流形計算方法計算出第一條Foliation弧長為ARC的流形軌道L1.

以L1為當前的參考軌道，順次選取初始圓上的下一點與x0相連構成初始流形段，進行一維流形軌道計算，使Foliation弧長為ARC.然后計算當前軌道上離散點與參考軌道上處于相同Foliation弧長處離散點的距離，若其中某個點對的距離大于網格最大半徑SIZEmax，則說明兩條軌道間的距離過大，需要在這兩條軌道之間再插入一條新的軌道.插入新軌道的方法為：在初始圓上兩條軌道的出發點的中點處插入一個網格點，然后以x0和該網格點的連線為初始流形段，計算一條Foliation弧長為ARC的軌道.比較新軌道與參考軌道的距離，若距離仍過大，則按上述方法繼續插入新的軌道，否則以新軌道為新的參考軌道，選取初始圓上的下一點進行一維流形軌道計算.

上述過程總的目的是計算出足夠多的一維流形軌道來對二維流形面進行覆蓋，所以不可避免地會產生冗余，有必要刪除一些相互之間距離過小的軌道.仿真過程中發現，一般情況下，相鄰軌道最末點之間的距離是兩條軌道之間最大的，所以我們以軌道末點之間的距離來代表軌道之間的距離.對于相鄰的三條軌道Li、Li＋1、Li＋2，若Li與Li＋1的距離小于網格最小半徑SIZEmin，并且Li與Li＋2的距離小于網格最大半徑SIZEmax，則刪除軌道Li＋1.

接下來就是對計算結果進行可視化.在每條軌道上，抽取Foliation弧長為k×step（k＝1，2，…）處的網格點構成新軌道，其中step規定了網格的大小，可以根據需要進行調整.由于原軌道上的網格點的Foliation弧長不一定恰好滿足step的整數倍條件，所以抽取過程中要根據原軌道上網格點的Foliation弧長進行新網格點的線性插值.抽取結束后，依次將各條新軌道上k＝1，2，…，的網格點用線連起來，以流形環的形式表示二維流形，所以同一環上的網格點的Foliation弧長是相同的.可以進一步進行三角形有限元網格化，將二維流形表示為一個曲面，相鄰環間的網格化如圖3所示，由于相鄰環上的網格點是一一對應的，所以操作起來也很方便.

二維穩定流形的計算除了一維流形軌道的計算方法不同，其余操作類似.

圖3 流形環之間的三角形網格化Fig.3 Triangulation between two neighboring circles

4 算例分析

仿真硬件環境為AMD phenom三核，2G內存的計算機，使用MATLAB編程.

4.1 三維Hénon映射的二維穩定流形

三維Hénon映射［20］的表達式

當M1＝1.4、M2＝0.2、B＝0.1時，該系統具有混沌吸引子，見圖4，這點與著名的二維Hénon映射類似.此時系統的不動點為x0＝（x，y，z），其中x ＝y ＝z＝0.883 9．F 在x0處的Jacobian矩陣為

矩陣A的3個特征值的模值分別為：λ1＝3.410 6、λ2＝0.038 6、λ3＝0.075 9，所以x0是一個雙曲不動點，且具有一維不穩定流形和二維穩定流形.

計算過程中采用的精度參數如下：δ＝0.1、βmax＝0.15、αmax＝0.15、αmin＝0.05、Δmax＝0.2、Δmin＝0.05、S＝1.4、SIZEmax＝1、SIZEmin＝0.5、step＝0.5．本文采用流形均勻增長的策略，而對于實際流形可能曲率變化不一致，增長步長越小，越能反映局部快速變化的細節，但計算量隨之增加，兩者是不可調和的矛盾，本文選取的增長步長是兼顧計算速度和計算精度的中間值.

為了說明本文提出的Foliation控制方法能夠實現二維流形在各個方向上的均勻增長，先計算一段ARC＝5的二維穩定流形，結果如圖5所示.其中不同灰度的環帶之間的距離為1，實線為實際計算過程中得到的一維流形軌道，軌道上的點則是抽取的Foliation弧長為step整數倍處的點，容易看出推廣的Foliation條件能夠很好地保證二維流形在各個方向上的均勻增長.計算圖5中的二維流形僅耗時1.25 s.

圖4 三維Hénon映射的混沌吸引子Fig.4 Chaotic attractor of 3D Hénon map

圖5 三維Hénon映射二維穩定流形的起始部分Fig.5 Localmanifold of a 3D Hénon map

下面計算ARC＝30的二維穩定流形，結果如圖6所示，耗時為108 s.其中不同灰度的環帶之間的距離為1.二維穩定流形的形狀類似于將一張圓形紙張進行了兩次對折.另外，由圖可見，當流形的弧長增加時，本文的算法仍然能夠很好地控制二維流形在各個方向上的增長速度.

圖6 三維Hénon映射的二維穩定流形Fig.6 Stablemanifold of a 3D Hénonmap（ARC＝30）

4.2 Lorenz系統的二維穩定流形

Lorenz系統是一個大氣對流模型，以其著名的蝶形混沌吸引子而備受關注.該模型的表達式

當σ＝10、ρ＝28、β＝8／3時，系統的蝶形吸引子如圖7（a）所示.該模型是用微分方程表示的，是一個連續系統，我們以差分的形式將其離散化

為了盡可能地保持原系統的性質，需要對采樣間隔T進行優化選取．通過實驗，發現當T＝0.01時，離散系統具有與原連續系統類似的混沌吸引子（見圖7（b）），而且系統的演化速度適中，比較好地保持了原系統的性質.

圖7 （a）原Lorenz系統混沌吸引子和（b）T＝0.01時，離散化后的Lorenz系統混沌吸引子Fig.7 （a）Chaotic attractor of a continuous Lorenz system；（b）Chaotic attractor of a discrete Lorenz system；

原點是離散化后的Lorenz系統的一個不動點，系統在原點處的Jacobian矩陣為

矩陣A的3個特征值的模值分別為：λ1＝0.771 7、λ2＝1.118 3、λ3＝0.973 3，所以x0是一個雙曲不動點，且具有一維不穩定流形和二維穩定流形，這點與原Lorenz系統類似.

計算過程中采用的精度參數如下：ARC ＝80、δ＝0.1、βmax＝0.15、αmax＝0.1、αmin＝0.05、Δmax＝0.1、Δmin＝0.001、S＝1.4、SIZEmax＝1、SIZEmin＝0.5、step＝0.8.

計算結果如圖8和圖9所示，計算耗時約為22分鐘.由圖易見本文的算法很好地控制了二維穩定流形在各個方向上的增長速度.該二維穩定流形具有很復雜的幾何結構，這也是Lorenz系統復雜動力學性能的一種反映.其中圖8以流形環的形式來表示二維穩定流形，不同灰度帶之間的距離為5，圖9以曲面的形式來表示二維穩定流形，不同灰度帶之間的距離為10.

圖8 以流形環來表示二維穩定流形Fig.8 Manifold covered by circles

圖9 以流形曲面來表示二維穩定流形Fig.9 Manifold covered by circles

圖10為系統的一維不穩定流形和二維不穩定流形.其中一維不穩定流形用曲線表示，弧長為1 000，計算耗時為20秒.從圖中可以觀察到，穩定和不穩定流形除了在原點處相交外，之后二維穩定流形面與蝶形的一維不穩定流形相互纏繞、卷曲，但是始終未相交，這也從一個側面反映了本文算法的準確性.

5 總結

研究離散動力系統的二維流形計算問題.首先從不變流形沿軌道具有導數傳遞這個特殊性質出發，提出了一種新的一維流形計算方法，計算過程中以一維流形軌道的切線為參考基準來檢查精度，較文獻［10］中的控制條件更為簡潔，同時能夠保證計算的精確性.更為重要的是，本文的算法可以同時用于一維穩定流形和一維不穩定流形的計算.以一維流形計算方法為基礎，提出了二維流形計算方法.針對文獻［17］中提出的Foliation條件只能應用于流形面的局限，通過對葉形重新進行定義，提出了基于流形軌道的Foliation條件，仿真結果表明，推廣的Foliation條件能夠很好地控制二維流形在各個方向上的增長速度.文獻［17］的算法由于默認新網格點的原像處于已有流形的邊界上，所以應用中會有所局限，而本文的方法是以一維流形軌道曲線的計算為基礎的，應用中不存在這個限制.相對于文獻［18］中的算法，本文顯然在均衡二維流形在各個方向的增長速度上更具優勢，另外，文獻［18］中的算法僅是針對于二維不穩定流形的計算提出的，而本文的算法可以同時計算二維穩定流形和二維不穩定流形，文中的仿真例子說明了這一點，所以本文的算法在適用范圍上比文獻［18］中的算法更廣.

由于本文算法在進行插值時直接從初始環處開始計算軌道，這樣就避免了插值點原像的搜索問題，而且插值點之后的軌道也很容易計算，這樣做的一個冗余之處在于，會在二維流形上Foliation弧長較短處產生大量的離散點，這也是本文算法以后改進的一個方向.但是值得強調的一點是，對離散動力系統的二維流形進行網格點插值時，需要搜索插值點原像的位置，而這種搜索是在二維流形曲面上進行的，所以操作起來較直接計算軌道更為復雜，這也是本文采取直接計算軌道法的原因.使用本文的方法計算同宿軌道和異宿軌道的前提是必須已知同宿橫截點或者異宿橫截點以及其對應的增長方向，這是計算同宿軌道和異宿軌道比較困難的問題.

圖10 離散化后Lorenz系統的一維不穩定流形和二維穩定流形Fig.10 Stable and unstablemanifold of discrete Lorenz system

［1］Johnson ME，Jolly MS，Kevrekidis IG.Two-dimensional invariantmanifolds and global bifurcations：Some approximation and visualization studies［J］.Numer Algorithms，1997，14：125－140.

［2］Guckenheimer J，Worfolk P，Dynamical systems：Some computational problems［M］.Kluwer Academic Publishers.1993：241－277.

［3］Henderson ME，Siam J.Computing invariantmanifolds by integrating fattened trajectories［J］.Appl Dyn Syst，2005，4：832－882.

［4］Krauskopf B，Osinga H M.Two-dimensional globalmanifolds of vector fields［J］.Chaos，1999，9：768－774.

［5］Krauskopf B，Osinga H M，Doedel E J，Henderson ME，Guckenheimer J，Vladimirsky A，Dellnitz M，Junge O.A survey of methods for computing（un）stablemanifolds of vector fields［J］.Bifur Chaos Appl Sci Engrg，2005，15：763－791.

［6］Guckenheimer J，Vladimirsky A，Siam J.A fastmethod for approximating invariantmanifolds［J］.Appl Dyn Syst，2004，3：232－260.

［7］Li Q D，Yang X S.A new algorithm for computation of two-dimensional unstablemanifolds and its applications［J］.Acta Phys Sin，2010，59（3）：1416－1422.

［8］You Z，Kostelich EJ，Yorke JA，Calculating stable and unstablemanifolds［J］.Int JBifurc Chaos Appl Sci Engrg，1991，1：605－623.

［9］SimóC.On the analytical and numerical approximation of invariantmanifolds［M］∥Benest D，FroeschléC.Les Méthodes modernes de lamécanique céleste goutelas，1989，285－330.

［10］Parker TS，Chua LO.Practical numerical algorithms for chaotic systems［M］.Berlin：Springer，1989：131－140.

［11］Hobson D.An efficientmethod for computing invariantmanifolds［J］.JComput Phys，1991，104：14－22.

［12］Krauskopf B，Osinga H M.Growing unstable manifolds of planar maps［EB／OL］.1997，http：／／www.ima.umn.edu／preprints／OCT97／1517.ps.gz.

［13］Krauskopf B，Osinga H M.Growing 1D and quasi-2D unstablemanifolds ofmaps［J］.JComput Phys，1998，146（1）：406－419.

［14］England JP，Krauskopf B，Osinga H M.Computing one-dimensional stablemanifolds and stable sets of planarmaps without the inverse［J］.Siam JAppl Dyn Syst，2004，3（2）：161－190.

［15］Dellnitz M，Hobmann A.A subdivision algorithm for the computation of unstablemanifolds and global attractors［J］.Numer Math，1997，75：293－317.

［16］Fundinger D.Toward the calculation of higher-dimensional stable manifolds and stable sets for noninvertible and piecewisesmooth maps［J］.JNonlinear Sci，2008，18：391－413.

［17］Krauskopf B，Osinga H M.Globalizing two-dimensional unstablemanifolds ofmaps［J］.Int JBifurc Chaos Appl Sci Engrg，1998，8：483－503.

［18］Li Q D，Zhou L，Zhou HW.Computation for two-dimension unstablemanifold ofmap［J］.Journal of Chongqing University of Posts and Telecommunications（Natural Science Edition），2010，22（3）：339－345.

［19］Palis J，MeloWD.Geometric theory of dynamical systems［M］.New York：Springer-Verlag，1982.

［20］Gonchenko SV，Ovsyannikov II，Simo C，Turaev D.Three-dimensional Hénon-likemaps and wild Lorenz-like attractors［J］. Internat JBifur Chaos Appl Sci Engrg，2005，15（11）：3493－3508.

［21］Li Q D，Yang X S.Computation of two dimensional unstablemanifold［J］.Chinese Journal of Computational Physics，2005，22（6）：549－554.

［22］Timothy JH，Li Q D，Cheng R B.Wrinkling behavior of highly stretched rectangular elastic films via parametric global bifurcation［J］.Journal of Nonlinear Science，2013，23：777－805.

［23］Sun H，Fan Y，Li H，et al.Computation of two-dimensional invariant manifolds with radial growth factor［J］.Chinese J Comput Phys，2011，28（4）：621－625.

［24］Li H，Fan Y，Zhang J.Computing stable and unstablemanifolds of typical chaoticmaps［J］.Chinese JComput Phys，2011，28（6）：927－932.

Grow ing 2D Manifold of Discrete Dynam ical System Based on Foliation Condition

JIA Meng
（Department of Electrical Engineering，Xinxiang College，Xinxiang 453003，China）

An algorithm for computing 2D stable and unstablemanifolds of hyperbolic fixed points of discrete dynamical systems is shown.With the fact that Jacobian transports derivative along orbitof an invariantmanifold，an algorithm for computing 1Dmanifold is proposed.Themesh point is located with a Prediction-Correction schemewhich reduces searching time and at the same time gives rise to a simplified accuracy condition.2D manifold is computed by covering it with orbits of 1D sub-manifold.A generalized Foliation condition is used to guarantee that2D manifold is growing equally along orbits of 1D sub-manifold in different directions.Performance of the algorithm is demonstrated with hyper chaotic 3D Hénonmap and Lorenz system.

discrete dynamical system；stable manifold；unstable manifold；derivative transportation；3D Hénon map；Lorenz system；chaotic attractor

date：2013－07－22；Revised date：2013－11－26

O59

A

1001-246X（2014）04-0495-10

2013－07－22；

2013－11－26

河南省重大科技攻關項目（112102210014）資助

賈蒙（1981－），男，河南新鄉，博士，從事非線性流形計算和非線性信號處理研究，E-mail：tianshi_cd＠163.com