賈 蒙

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

動力系統(tǒng)按照其狀態(tài)變量隨時間變化的連續(xù)性可分為兩類：連續(xù)動力系統(tǒng)和離散映射動力系統(tǒng)。其中連續(xù)動力系統(tǒng)一般可用微分方程的形式來表達(dá),而離散動力系統(tǒng)則一般用點到點的映射來表示,這也是系統(tǒng)演化的一種表達(dá)方式,可看作是對連續(xù)動力系統(tǒng)采樣后的結(jié)果?？茖W(xué)研究及生活實際中的很多過程都可以用一定的微分方程或者是離散映射來描述,本文側(cè)重于研究離散映射動力系統(tǒng)的計算。

本文主要討論離散映射動力系統(tǒng)中鞍型不動點的一維不變流形計算問題。穩(wěn)定流形和不穩(wěn)定流形在分析系統(tǒng)的動力學(xué)特性中起著非常重要的作用,它們充當(dāng)不同吸引子吸引域的邊界,將全空間劃分為多個具有不同動力特性的不變子空間,而當(dāng)穩(wěn)定與不穩(wěn)定流形相交時,就會引起同宿、異宿以及混沌等復(fù)雜動力學(xué)行為的出現(xiàn)。不變流形計算對于上述問題的研究都有著積極地促進(jìn)作用。

流形計算作為動力系統(tǒng)分析的一種幾何方法,一直以來都受到研究人員的重視。連續(xù)動力系統(tǒng)一般以微分方程的形式來表示,運用現(xiàn)代數(shù)值積分方法對相應(yīng)的初值直接進(jìn)行積分就能得到其一維流形,計算非常簡單,計算的精度也可以保證；其二維流形的計算雖然較一維流形的計算有了很大的挑戰(zhàn),但依舊出現(xiàn)了很多性能出眾的算法[1-3]。離散映射動力系統(tǒng)的流形計算相比較就沒有那么幸運了,即使是最為簡單的一維流形的計算也是困難重重,最常見的方法是先對系統(tǒng)在不動點的鄰域內(nèi)進(jìn)行線性化,得到局部流形,然后以此局部流形作為基域進(jìn)行迭代。這種算法操作簡單,而且計算速度快,但由于只能對基域上的有限個離散點進(jìn)行迭代,迭代過程中點的分離可能導(dǎo)致在基域上相距很近的兩個點經(jīng)迭代后得到的像之間距離變得很大,從而錯過流形的細(xì)節(jié),使得計算的精度變差。You等[4]，Simó[5]，Parker等[6],Hobson[7]給出了基域迭代的改進(jìn)算法。Krauskopf等[8-9]給出了一種有別于基域迭代的新算法,他們一步步來增長流形,并結(jié)合了曲率控制技術(shù),流形的增長速度由流形的局部曲率來決定,該算法能夠有效展示流形的細(xì)節(jié),與此思想類似,他們在文獻(xiàn)[10]中提出的算法能夠在逆映射不能顯式表達(dá)的情況下計算系統(tǒng)的穩(wěn)定流形,而之前提到的大部分算法都是把系統(tǒng)的穩(wěn)定流形當(dāng)作其逆映射的不穩(wěn)定流形來計算的。另外還有Dellnitz等[11]提出的細(xì)分法,該算法先用“盒子”覆蓋所要研究的區(qū)域,然后逐步對盒子進(jìn)行細(xì)分,從而得出所要計算的流形；Fundinger[12]也提出了與此思想類似的一種算法。細(xì)分法同時還可用于離散動力系統(tǒng)二維流形的計算,Krauskopf等[13]提出的算法也可用于離散動力系統(tǒng)二維流形的計算。

本文提出的算法也是從基域迭代這個思想出發(fā)的,同時借鑒了參考文獻(xiàn)[7,10]中曲率控制及單次映射的長處,該算法易于編程實現(xiàn),同時也可以并行進(jìn)行計算。

1 離散映射系統(tǒng)流形計算的問題分析

對于系統(tǒng)(1),

xn+1=F(xn)

(1)

假設(shè)其存在不動點x0,滿足F(x0)=x0,要研究系統(tǒng)在x0附近的特性,一般都要先對系統(tǒng)進(jìn)行線性化,將式(1)轉(zhuǎn)化為線性映射x→Ax,其中A為x0處的雅各比矩陣A=DF(x0)=[?fi/?xj](x0)。若矩陣A的特征值的模都不等于1,那么x0就是一個雙曲不動點；從幾何角度來說,就是在復(fù)平面上,矩陣A的特征值均處在單位圓內(nèi)或者單位圓外,而沒有落在單位圓上的。其中處于單位圓內(nèi)的特征值叫做穩(wěn)定特征值,其對應(yīng)的特征向量{v1,v2,…,vl}張成穩(wěn)定特征空間Es；處于單位圓外的特征值叫做不穩(wěn)定特征值,它們對應(yīng)的特征向量{vl+1,vl+2,…,vn}張成不穩(wěn)定特征空間Eu。全空間En=Es⊕Eu。

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

(2)

(3)

全局穩(wěn)定和不穩(wěn)定流形定義為

(4)

(5)

現(xiàn)在的問題是,在不動點x0的鄰域U內(nèi),原映射系統(tǒng)F是否與其線性化x→Ax具有相同的性質(zhì)?？梢宰C明[15-16],對于雙曲不動點x0,非線性系統(tǒng)F與其線性化系統(tǒng)在x0附近是拓?fù)涔曹椀?從而具有相同的拓?fù)浣Y(jié)構(gòu)。

(6)

(7)

并且Wi∩Wj=?(i≠j),則由式(6)可知

(8)

式(8)表明,一維流形的計算可以分解為局部流形各個子區(qū)間的迭代值的并集,而各個子區(qū)間的計算互不影響,也就是說,可以將流形計算這個大任務(wù)劃分為m個小任務(wù)來執(zhí)行,這樣就使得并行計算的實現(xiàn)成為可能。

直接使用式(8)時,會使人產(chǎn)生疑問：各個子區(qū)間之間經(jīng)過迭代之后是否會有交集呢,如果有,則將造成重復(fù)計算,增加了不必要的計算量,那么如何進(jìn)行處理？

證明：首先定義一維不穩(wěn)定流形Wu(x0)上兩個點x,y之間的距離為

(9)

點與集合的距離定義為

(10)

(11)

其中:dg(W1,x0)=0距離最小，dg(Wm0,x0)距離最大。

由xi的任意性可知Fi(W1)∩Fj(Wk)≠?,(k=1,2,…,m0),類似可以證明在某些條件下,Fi(Wk1)∩Fj(Wk2)≠?,其中k1,k2=1,2,…,m0。

所以直接對鄰域U內(nèi)的局部流形進(jìn)行分割必將導(dǎo)致某些區(qū)間的像的重疊,證畢。

根據(jù)上文的討論,有必要對區(qū)間的確定進(jìn)行優(yōu)化,在保證劃分后的區(qū)間經(jīng)迭代產(chǎn)生的像沒有重疊(后邊為了編程的方便,相鄰區(qū)間的邊界是重合的)的同時,也要保證不能丟掉部分流形信息。這個問題看似非常復(fù)雜,但解決起來卻很容易,早期的基域迭代方法[4-6]就給出了一種非常聰明的辦法,雖然并沒有引起人們多大的注意。

他們的方法是：在Eu上沿不穩(wěn)定特征向量方向取距離x0為δ的一點p0,計算p0的像p1=F(p0),當(dāng)δ取較小值時，p0和p1之間的流形可近似用線段代替,以線段p0p1為迭代的初始段。下面將證明點p0、p1所在的流形分支上的所有點都可以通過F映射到線段p0p1上。

定理2對于離散映射函數(shù)(1),假設(shè)xk1和xk2處在雙曲平衡點x0一維不穩(wěn)定流形Wu(x0)的同一分支上的任意兩點,且滿足

dg(xk1,x0)>dg(xk2,x0)

則dg(Fk(xk1),x0)>dg(Fk(xk2),x0),其中k為任意整數(shù)。

證明：用反證法進(jìn)行證明

集合γFk(x),x={Wu(x0)上連續(xù)Fk(x)和x的曲線的所有點}。

如圖1所示,假設(shè)dg(Fk(xk1),x0)

則當(dāng)n→∞時,區(qū)間套的長度將逐漸減小到零,從而使得Fk(xkn)=xkn,也就是說Wu(x0)上除x0外,還存在不動點xkn,這與不穩(wěn)定流形Wu(x0)的定義相矛盾,所以假設(shè)不成立,原命題成立。證畢。

圖1 一維不穩(wěn)定流形某一分支上的映射

應(yīng)用定理2,可以很容易說明前面提到的問題。對于和點p0、p1處于同一流形分支的Wu(x0)上的點x,若dg(x,x0)>dg(p1,x0),則存在正整數(shù)k1使得

也就是說點F-k1(x)和F-(k1-1)(x)處在點p1的異側(cè),則由定理2可知

dg(F-k1(x),x0)=dg[F-1(F-(k1-1)(x)),x0]>

dg(F-1(p1),x0)=dg(p0,x0)

即F-k1(x)處于Wu(x0)上點p0、p1之間。同樣,當(dāng)dg(x,x0)

證明：令p0p1=W1∪W2∪…∪Wm0,且Wi∩Wj=?(i≠j),任取ain∈Wi、ajk∈Wj(i≠j),假定滿足dg(ain,x0)>dg(ajk,x0),則由定理2知dg(Fk(ain),x0)>dg(Fk(ajk),x0),由ain和ajk的任意性可知

所以Fk(Wi)∩Fk(Wj)=?。同理可證對Fk(p0p1)進(jìn)行任意分割后個子區(qū)間的像也不會重疊。由于任意一點都在p0p1上有原像,則由式(9)知各區(qū)間的像的并集就構(gòu)成了Wu(x0)。證畢。

另外,根據(jù)式(2)的性質(zhì),如果已知Fk(p0p1),則在計算Fk+1(p0p1)時,可直接將Fk(p0p1)看作當(dāng)前的迭代映射段進(jìn)行一次映射即可,而不需要再從p0p1開始進(jìn)行迭代,這樣做有利于減小計算量。

2 算法描述

假定當(dāng)前的迭代區(qū)域由離散點Mn={xn,1,xn,2,…,xn,k},其中下標(biāo)n表示第n次計算時的迭代域，k為迭代域中各離散點的順序號，k越大，dg(xnk,x0)也越大。應(yīng)當(dāng)注意到,當(dāng)前迭代區(qū)域的首末兩點滿足xn,k=F(xn,1),所以我們從Mn的第二個點開始計算映射值,但為了保證前后段之間的映射對應(yīng)關(guān)系,仍將xn,k作為Mn+1中的第一個點。由于迭代過程中可能導(dǎo)致點的分離,使得Mn上兩個相鄰點的經(jīng)映射后得到的像的距離變得很大,從而不能很好地逼近實際流形,為了保證計算的精度,在計算完Mn中所有離散點的映射點之后,要對映射點進(jìn)行精度檢查。采用文獻(xiàn)[9-10]中提出的曲率約束條件。

考慮軌道上相鄰三個離散點之間的角度,如圖2所示。

由圖得

(12)

其中

(13)

圖2 曲率約束條件

(14)

現(xiàn)在檢查α是否滿足下列約束條件

(15)

利用上述曲率約束條件,我們可以確定Mn上兩個相鄰點的經(jīng)F映射后產(chǎn)生的新點是否滿足精度條件。假定已經(jīng)得到了Mn中第i+1個點xn,i+1的映射點F(xn,i+1),檢驗其精度時發(fā)現(xiàn)其不滿足條件,則先保存F(xn,i+1)的坐標(biāo)值,然后需要在F(xn,i+1)之前插入其他離散點,直至滿足精度條件。需要說明的是，xn+1,j和xn+1,j+1之間可能需要插入不止一個點才能滿足精度條件,而且所插入的點也要檢查其精度條件。由于Mn和Mn+1之間具有映射對應(yīng)關(guān)系,所以當(dāng)需要在F(xn,i)和F(xn,i+1)之間插入其他離散點時,其原像都近似處于點xn,j和xn,j+1的連線上,因此無需再搜索插入點原像的位置。關(guān)于上述映射點的計算以及點的插入算法我們將在后面詳細(xì)介紹。

新的迭代段：

Mn+1={F(xn,1),F(xn,1,1),…,

F(xn,1,a),F(xn,2),…,F(xn,i,j),…,

F(xn,k-1),F(xn,k-1,),…,F(xn,k-1,b),F(xn,k)}

其中:F(xn,k)為Mn中各離散點的映射值，F(xiàn)(xn,i,j)為在檢查精度過程中在F(xn,i)和F(xn,i+1)之間插入的點。

這時候我們對Mn+1中的點再進(jìn)行一次精度檢查：若αmin<α、 (Δα)min<Δα并且Δ<Δmin,則說明離散點的間距過小,就要刪除一個點；當(dāng)αmin<α、 (Δα)min<Δα并且Δ>Δmin時,認(rèn)為離散點的間距過大,就使用線性插值插入一個點。為了保證Mn+1上首尾兩點間的映射對應(yīng)關(guān)系,即使首末兩點滿足了刪除條件,依舊保留,而不進(jìn)行刪除(可以證明,刪除后將導(dǎo)致流形部分信息丟失)。本文在控制精度時綜合運用了曲率約束和距離控制兩種技術(shù),這樣利于保證在流形變化比較緩慢的區(qū)域相鄰離散點間的距離不會太大,這是與文獻(xiàn)[7-8]的不同之處。

精度檢查結(jié)束后,將Mn+1加入到離散序列點集合M中(由于Mn+1中的第一個點與M的最末點相同,所以加入時要進(jìn)行相應(yīng)的覆蓋操作)。判斷弧長arcl是否達(dá)到總弧長ABC,否的話把Mn+1作為新的當(dāng)前迭代段繼續(xù)進(jìn)行計算,是則結(jié)束。

3 點的插入算法及并行計算的詳細(xì)討論

先討論點的插入算法。依次計算Mn中各離散點的映射值,并檢查其精度條件,不滿足精度條件時就需要插入其它離散點。假設(shè)現(xiàn)在要在點F(xn,i)和F(xn,i+1)間進(jìn)行點的插入操作。取區(qū)間xn,ixn,i+1的中點xn,i+1/2,然后計算其映射點F(xn,i+1/2),并判斷該點的位置是否滿足精度要求：若不滿足,則保存F(xn,i+1/2)的值,再取區(qū)間xn,ixn,i+1/2的中點,計算相應(yīng)的映射點,檢查精度,就這樣循環(huán)下去直至滿足條件；若F(xn,i+1/2)的位置滿足精度條件,則將F(xn,i+1/2)加入到Mn+1的序列中,然后檢查F(xn,i+1)的精度條件,若不滿足,與前邊類似,要找區(qū)間xn,i+1/2xn,i+1的中點,計算相應(yīng)的映射點,檢查精度條件,就這樣一直折半映射下去直至滿足精度條件,并將插入的離散點加入到Mn+1的序列中。需要注意的是,每次區(qū)間分半后都牽扯到兩個小區(qū)間上的精度檢查。

上面的敘述看起來非常繁瑣,但歸納起來只有三個步驟：① 計算區(qū)間xn,ixn,i+1中點xn,i+1/2的映射值F(xn,i+1/2)并保存；② 檢查區(qū)間的前半段F(xn,i)F(xn,i+1/2)上的精度,滿足則結(jié)束,不滿足則在F(xn,i)F(xn,i+1/2)用相同方法插入足夠點直至滿足精度條件；③ 檢查區(qū)間的后半段F(xn,i+1/2)F(xn,i+1)上的精度滿足則結(jié)束,不滿足則在F(xn,i+1/2)F(xn,i+1)用相同方法插入足夠點直至滿足精度條件。

由于所有的點的插入都遵循了“先找區(qū)間中點并計算映射點然后檢查精度”同樣的過程,基本思想是在兩個不滿足精度條件的離散點中間插入一個離散點使其滿足精度,所以上面的三個步驟可用遞歸程序來實現(xiàn),形式簡潔而且編程實現(xiàn)也容易。

4 仿真實例

以下兩個例子都是在Matlab環(huán)境下編程實現(xiàn),計算機配置為AMD三核2.3 GHz、2 G內(nèi)存。為了便于比較計算的結(jié)果,先介紹文獻(xiàn)[8]中的一個例子。

4.1 shear map

計算時采用的參數(shù)為λu=2、λs=0.4、δ=0.001，αmax=0.3、αmin=0.2、 Δmax=0.01、 Δmin=0.000 1、 (Δα)max=10-5、 (Δα)min=10-6。計算結(jié)果如圖3所示。耗時對比結(jié)果示于表1圖3中(a)耗時12 s,(b)耗時3 s,(c)耗時16 s,(d)耗時約88 s。圖中只畫出了橫軸正半軸范圍內(nèi)的不穩(wěn)定流形的分支,另一分支處于負(fù)橫軸上,與理論分析一致。

表1 Shear map不同算法耗時對比結(jié)果

圖3 Shear map在原點處的一維不穩(wěn)定流形

4.2 GHM (廣義Hénon映射)

廣義Hénon映射[17-18]經(jīng)常出現(xiàn)在同宿分岔的理論分析中,其表達(dá)式為

F2在x0處的線性化矩陣為：

其中:a11=-b+Ry，a12=-2y+Rx；

a21=-2(a-bx-y2+Rxy)×

(-bx+Ry)+Ry(-bx+Ry);

a22=-b-2(a-bx-y2+Rxy)×

(-2y+Ry)+Ra-Rbx-3Ry2+2

A在x0處具有特征值λs=0.245 588 874 766 114、λu=2.762 615 644 582 639,所以F2在x0處具有一維不穩(wěn)定流形Wu(x0)，λu對應(yīng)的特征向量為Vu=(0.515 5,0.856 9)T。

計算過程中參數(shù)設(shè)置如下：Arc=10，δ=0.001，αmax=0.3、αmin=0.2、 Δmax=0.2、 Δmin=0.000 1、 (Δα)max=10-3、 (Δα)min=10-4。計算結(jié)果如圖4所示。圖中方形標(biāo)注的是不動點x0。

圖4 廣義Hénon映射的一維不穩(wěn)定流形

圖5 仿真計算結(jié)果

圖5(a)是文獻(xiàn)[17]計算得出的結(jié)果,圖5(b)中實線為本文中方法的計算結(jié)果,由圖可觀察到,兩種方法的計算結(jié)果比較吻合。值得一提的是,圖5(a)計算耗時12.5 s,而圖4耗時0.09 s,相比較而言,后者計算速度約為前者的140倍。

5 精度討論

(16)

舍掉高階無窮小量,誤差:

E=F(xu+εi)-F(xu)=F′(xu)εi=Aεi

(17)

對于二維相空間

(18)

(19)

以A的兩個特征向量u,v為基張成二維空間R2,

(20)

對于一維不穩(wěn)定流形Wu(x0),初始誤差ε0主要處于方向v上,而在u上的分量很小,所以ε0≈m2v,則經(jīng)過一次映射之后,誤差

ε1=Aε0≈Am2v=m2λ2v

(21)

所以初始誤差經(jīng)迭代后將變小。對于任意一點xu處的插值誤差εi,我們無法直接寫出類似于式(18)的表達(dá)式,所以通過式(20)來分析：

可見插值誤差εi經(jīng)映射后造成的傳播誤差依然保持有界。如果xu處的線性化矩陣具有系統(tǒng)不動點x0處線性化矩陣A類似的性質(zhì),則可以說明傳播誤差也是減小的,這點還有待證明。

可以確定的是,當(dāng)提高計算時所采用的精度參數(shù)時，ε0和εi都會隨之減小,但是嚴(yán)格的全局誤差界還是比較難以確定。一種檢驗辦法是,不斷提高計算的精度參數(shù),然后比較各次的計算結(jié)果,當(dāng)解的變化小于某一精度值時,計算結(jié)束。

下面討論第4節(jié)中例子的精度,改變參數(shù)δ,然后比較兩次計算所得流形上處于相等弧長處的點的距離,見表2。

表2 Shear map初始步長δ取不同值時的計算結(jié)果(c=0.5)

由表2可知兩次計算的結(jié)果近似相同,偏差較小,并且隨著弧長的增加,誤差仍舊保持在原數(shù)量級,而沒有顯著增加。

為了比較串并計算的結(jié)果,在并行結(jié)果上選取弧長為0.25n處的點,然后計算這些點串行結(jié)果中流形的距離。如圖6所示。

圖6 誤差分析

由圖6(a)可見,并行結(jié)果上的離散點都分布在串行結(jié)果上,由圖6(b)可知離散點與流形線的偏差始終保持在10-6,說明串行計算和并行計算結(jié)果的差異是很小的。

6 結(jié) 論

本文提出了一種計算離散動力系統(tǒng)鞍點的一維不穩(wěn)定流形的并行快速算法,其基本思想源于基域迭代,并結(jié)合了文獻(xiàn)[7-8]中在曲率控制上計算速度快的優(yōu)點；此外,本文對局部流形進(jìn)行了剖分,說明了如何選取局部流形,從而保證剖分后的區(qū)間經(jīng)映射之后的映像不會重疊,并且不會漏掉流形段；從理論上說明了并行計算的可能性,并討論了并行計算時所要考慮的一些細(xì)節(jié)問題。當(dāng)要求計算很大弧長的流形時,并行計算是很有優(yōu)勢的。

除了可以并行計算,新算法還突出編程的簡便性,文中介紹了當(dāng)離散點間不滿足精度條件時的一種遞歸算法,將插值問題簡化為在不滿足精度條件的兩個點之間插入一個點使其滿足精度,插值的目的明確,編程簡潔,易于實現(xiàn)。此外,由前面的例子可知,即使是串行計算,新方法在計算速度上也具有優(yōu)勢。

[1]Krauskopf B，Osinga H M，Doedel E J，et al.A survey of methods for computing (un)stable manifolds of vector fields[J].Bifur.Chaos.Appl.Sci.Engrg，2005，15(3):763-791.

[2]Guckenheimer J，Vladimirsky A.A fast method for approximating invariant manifolds[J].Appl.Dyn.Sys，2004，3(3): 232-260.

[3]Henderson M E.Computing invariant manifolds by integrating fattened trajectories[J].Appl.Dyn.Sys，2005，4(4):832-882.

[4]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.

[5]Sim C.On the analytical and numerical approximation of invariant manifolds.In: Benest，D.，F(xiàn)roeschl ,C.(eds.)Les M thodes Modernes de la M canique C leste，Goutelas，1989:285-330.

[6]Parker T S，Chua L O.Practical numerical algorithms for chaotic systems.Springer，Berlin,1989.

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

[8]Krauskopf B，Osinga H M，Growing unstable manifolds of planar maps，1517，1997，http://www.ima.umn.edu/preprints/OCT97/1517.ps.gz.

[9]Krauskopf B，Osinga H M.Growing 1D and quasi-2D unstable manifolds of maps[J].J.Comput.Phys,1998b,146: 406-419.

[10]England J P，Krauskopf B，Osinga H M.Computing one-dimensional stable manifolds and stable sets of planar maps without the inverse[M].SIAM J.Appl.Dyn.Syst,2004,2:161-190.

[11]Dellnitz M，Hohmann A.A subdivision algorithm for the computation of unstable manifolds and global attractors[J].Numer.Math,1997,75: 293-317.

[12]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.

[13]Krauskopf B，Osinga H M.Globalizing two-dimensional unstable manifolds of maps[J].Int.J.Bifurc.Chaos,1998a,8(3): 483-503.

[14]Palis J，Melo W D.Geometric theory of dynamical systems[M].Springer，1982.

[15]Berglund N.Geometrical theory of dynamical systems.http://arxiv.org/abs/math/0111177v1.

[16]Jia M，Yang Y，Tian W J.A fast computing method to distinguished hyperbolic trajectory for non-autonomous system[J].Chinese Physics B，2011，20(3): 034701-1-034701-5.

[17]Govaerts W，Khoshsiar G R，Kuznetsov Y A，et al.A toolbox for continuation and bifurcation of cycles of maps，http://www.matcont.UGent.be,2006.

[18]Ghaziani R K，Govaerts W，Yu A，et al.Numerical continuation of connecting orbits of maps in MATLAB[J].J.Diff.Eqns.Appl，2009,15(8):849-875.