映射動(dòng)力系統(tǒng)一維流形并行計(jì)算方法

賈 蒙

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

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

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

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

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

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

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

xn+1=F(xn)

(1)

假設(shè)其存在不動(dòng)點(diǎn)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就是一個(gè)雙曲不動(dòng)點(diǎn)；從幾何角度來說,就是在復(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的一個(gè)不動(dòng)點(diǎn),則在x0的鄰域U內(nèi)存在局部穩(wěn)定和不穩(wěn)定流形：

(2)

(3)

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

(4)

(5)

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

(6)

(7)

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

(8)

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

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

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

(9)

點(diǎn)與集合的距離定義為

(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ū)間的邊界是重合的)的同時(shí),也要保證不能丟掉部分流形信息。這個(gè)問題看似非常復(fù)雜,但解決起來卻很容易,早期的基域迭代方法[4-6]就給出了一種非常聰明的辦法,雖然并沒有引起人們多大的注意。

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

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

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的曲線的所有點(diǎn)}。

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

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

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

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

也就是說點(diǎn)F-k1(x)和F-(k1-1)(x)處在點(diǎn)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)上點(diǎn)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)行任意分割后個(gè)子區(qū)間的像也不會(huì)重疊。由于任意一點(diǎn)都在p0p1上有原像,則由式(9)知各區(qū)間的像的并集就構(gòu)成了Wu(x0)。證畢。

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

2 算法描述

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

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

由圖得

(12)

其中

(13)

圖2 曲率約束條件

(14)

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

(15)

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

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

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

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

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

上面的敘述看起來非常繁瑣,但歸納起來只有三個(gè)步驟：① 計(jì)算區(qū)間xn,ixn,i+1中點(diǎn)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)用相同方法插入足夠點(diǎn)直至滿足精度條件；③ 檢查區(qū)間的后半段F(xn,i+1/2)F(xn,i+1)上的精度滿足則結(jié)束,不滿足則在F(xn,i+1/2)F(xn,i+1)用相同方法插入足夠點(diǎn)直至滿足精度條件。

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

4 仿真實(shí)例

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

4.1 shear map

計(jì)算時(shí)采用的參數(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。計(jì)算結(jié)果如圖3所示。耗時(shí)對比結(jié)果示于表1圖3中(a)耗時(shí)12 s,(b)耗時(shí)3 s,(c)耗時(shí)16 s,(d)耗時(shí)約88 s。圖中只畫出了橫軸正半軸范圍內(nèi)的不穩(wěn)定流形的分支,另一分支處于負(fù)橫軸上,與理論分析一致。

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

圖3 Shear map在原點(diǎn)處的一維不穩(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。

計(jì)算過程中參數(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。計(jì)算結(jié)果如圖4所示。圖中方形標(biāo)注的是不動(dòng)點(diǎn)x0。

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

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

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

5 精度討論

(16)

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

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

(17)

對于二維相空間

(18)

(19)

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

(20)

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

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

(21)

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

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

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

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

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

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

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

圖6 誤差分析

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

6 結(jié) 論

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

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