高允報(bào)，孫玉泉

(北京航空航天大學(xué)數(shù)學(xué)科學(xué)學(xué)院，北京 100191)

Lyapunov指數(shù)(LE)可以判斷動(dòng)力系統(tǒng)的混沌狀態(tài)，在研究混沌運(yùn)動(dòng)特性中起著非常重要的作用，是衡量動(dòng)力學(xué)特性的一個(gè)重要定量指標(biāo)，它表征了系統(tǒng)在相空間中相鄰軌道間收斂或發(fā)散的平均指數(shù)率[1].系統(tǒng)的混沌狀態(tài)可以通過LE的正負(fù)來判斷，所以研究如何計(jì)算LE就非常必要.

近年來國(guó)內(nèi)外學(xué)者對(duì)LE的計(jì)算方法做了很多研究.嚴(yán)雯等[2]利用定義法求解了Logistic模型的最大Lyapunov指數(shù)(LLE)；Wolf等[3]利用Wolf法計(jì)算出幾種混沌系統(tǒng)的LLE；Rosenstein等[4]提出小數(shù)據(jù)量法來計(jì)算LLE.然而定義法只適用于計(jì)算已知系統(tǒng)的LE；Wolf法計(jì)算LLE時(shí)過程較為復(fù)雜且結(jié)果往往不精確；小數(shù)據(jù)量法在不知系統(tǒng)方程而僅可獲得離散數(shù)據(jù)的情況下計(jì)算LLE時(shí)，對(duì)于線性區(qū)域的判斷往往通過直觀觀察，導(dǎo)致結(jié)果誤差較大.本文將模糊C均值聚類算法用于小數(shù)據(jù)量法線性區(qū)域的選擇，提出了模糊C均值聚類小數(shù)據(jù)量法，提高了計(jì)算時(shí)間序列LLE的精度，且通過適當(dāng)增加每次迭代的離散時(shí)間步長(zhǎng)加快了計(jì)算LLE的速度.

1 預(yù)備知識(shí)

定義1.1Lyapunov指數(shù)：

(1)

經(jīng)過n次迭代后，這兩點(diǎn)的距離變?yōu)?/p>

δxn=|fn(x0+δx0)-fn(x0)|.

(2)

定義LE如下[5-6]：

(3)

對(duì)于一般的n維動(dòng)力系統(tǒng)，定義LE如下：設(shè)一個(gè)n維動(dòng)力系統(tǒng)xn+1=F(xn)(Rn→Rn)，將系統(tǒng)的初始條件取為一個(gè)無窮小的n維小球，半徑取為ε(0)，由于演化過程中的自然變形，小球?qū)⒆兂蓹E球，見圖1.將橢球上所有主軸按照其長(zhǎng)度排列，那么第i個(gè)LE根據(jù)第i個(gè)主軸的長(zhǎng)度εi(t)的增加速率[7-8]定義為

(4)

最大的λi稱為L(zhǎng)LE.系統(tǒng)的混沌狀態(tài)可以通過LLE的正負(fù)來判斷，當(dāng)LLE大于零時(shí)系統(tǒng)處于混沌狀態(tài)；當(dāng)LLE小于零時(shí)系統(tǒng)穩(wěn)定.

圖1 3維球面隨時(shí)間的演化

根據(jù)定義計(jì)算LE的方法稱為定義法，由于定義法計(jì)算LE時(shí)只適用于給定的系統(tǒng)，不能求時(shí)間序列[7-9]的LLE，因此Rosenstein提出了小數(shù)據(jù)量法.

定義1.2小數(shù)據(jù)量法：

對(duì)于時(shí)間序列{x1，x2，…，xN}的LLE可以用小數(shù)據(jù)量法，對(duì)于時(shí)間序列通過計(jì)算得[10]：設(shè)嵌入維數(shù)m，時(shí)間延遲為τ，平均周期為P，則重構(gòu)相空間為

Xi=[xi，xi+τ，…，xi+(m-1)τ]，i=1，2，…，M.

(5)

在重構(gòu)的相空間中尋找每個(gè)參考點(diǎn)Xj的最近鄰點(diǎn)Xj′，記

djj′(0)=min‖Xj-Xj′‖.

(6)

其中|j-j′|>P，P為時(shí)間序列的平均周期，可以通過基本軌道上每個(gè)點(diǎn)的平均發(fā)散速率計(jì)算.

對(duì)于每個(gè)參考點(diǎn)Xj，計(jì)算出其與最近鄰點(diǎn)Xj′的第i個(gè)離散時(shí)間步長(zhǎng)后的距離djj′(i)為

djj′(i)=min‖Xj+i-Xj′+i‖，i=1，2，…，min(M-j，M-j′).

(7)

假定參考點(diǎn)Xj與最近鄰點(diǎn)Xj′的指數(shù)發(fā)散率為λ1，那么

djj′(i)=Cjeλ1(i·Δt)，Cj=djj′(0)，

(8)

上式兩邊取對(duì)數(shù)得

lndjj′(i)=lnCj+λ1(i·Δt).

(9)

由(9)式可以看出lndjj′(i)與變量i滿足線性關(guān)系，其曲線斜率為λ1Δt.

因此固定i，對(duì)所有的lndjj′(i)求平均再除以Δt，得到平均發(fā)散程度指數(shù)

(10)

2 模糊C均值聚類小數(shù)據(jù)量法

2.1 模糊C均值聚類算法

在眾多模糊聚類算法中，模糊C-均值算法通過優(yōu)化目標(biāo)函數(shù)得到每個(gè)樣本點(diǎn)對(duì)所有類中心的隸屬度，從而決定樣本點(diǎn)的類屬以達(dá)到自動(dòng)對(duì)樣本數(shù)據(jù)進(jìn)行分類的目的.過程如下：

(11)

(12)

采用迭代的方法，通常情況下取初始隸屬度矩陣為每行為隨機(jī)單位向量的矩陣，b=2，最大迭代次數(shù)為100，隸屬度最小變化量為10-6，求解(11)和(12)式，直至滿足收斂條件，得到最優(yōu)解.

2.2 模糊C均值聚類算法選擇線性區(qū)域

在小數(shù)據(jù)量法計(jì)算LLE時(shí)將模糊C均值聚類算法用于選取線性區(qū)域，具體過程如下：

(1) 運(yùn)用模糊C均值聚類算法將y(i)數(shù)據(jù)分為不飽和區(qū)域(y(i)達(dá)到基本穩(wěn)定之前)和飽和區(qū)域(y(i)在一條水平直線附近做微小波動(dòng))兩類，如圖2所示.

(2) 對(duì)于不飽和區(qū)域需要選出線性區(qū)域，用二階差分ddy(i)近似替代二階導(dǎo)數(shù)，對(duì)于不飽和區(qū)域的二階差分(i，ddy(i))將其分成3類，如圖3所示.結(jié)合不飽和區(qū)域平均發(fā)散程度指數(shù)的圖像與不飽和區(qū)域的二階差分的圖像可以看出，當(dāng)i較小時(shí)離散程度較大，當(dāng)i較大時(shí)線性化程度較低.因此將不飽和區(qū)域的y(i)-i分成3類，如圖4所示，選取第2類作為最終選擇的線性區(qū)域，其斜率就是LLE.

圖2 平均發(fā)散程度指數(shù)

圖3 不飽和區(qū)域平均發(fā)散程度指數(shù)二階差分分類

2.3 模糊C均值聚類小數(shù)據(jù)量法

對(duì)于時(shí)間序列{y1，y2，…，yN}，嵌入維數(shù)my，時(shí)間延遲為τy，平均周期為Py，則重構(gòu)相空間為

Yi=[yi，xi+τy，…，yi+(my-1)τy]，i=1，2，…，My.

(13)

在重構(gòu)的相空間中尋找每個(gè)參考點(diǎn)Yj的最近鄰點(diǎn)Yj′，有

Djj′(0)=min‖Yj-Yj′‖.

(14)

其中|j-j′|>Py，Py為時(shí)間序列的平均周期，可以通過基本軌道上每個(gè)點(diǎn)的平均發(fā)散速率計(jì)算.

對(duì)于每個(gè)參考點(diǎn)Yj，計(jì)算出其與最近鄰點(diǎn)Yj′的第ni(n=1，2，3，…)個(gè)離散時(shí)間步長(zhǎng)后的距離為

Djj′(ni)=‖Yj+ni-Yj′+ni‖，

(15)

假定參考點(diǎn)Yj與最近鄰點(diǎn)Yj′的指數(shù)發(fā)散率為λ，那么

Djj′(ni)=Ajeλ(ni·Δt)，Aj=Djj′(0)，

(16)

上式兩邊取對(duì)數(shù)得

lnDjj′(ni)=lnAj+λ(ni·Δt).

(17)

由(17)式可以看出lnDjj′(ni)與i滿足線性關(guān)系，其曲線斜率為nλΔt.

因此固定i，對(duì)所有的lnDjj′(ni)求平均再除以Δt，得到平均發(fā)散程度指數(shù)

(18)

3 數(shù)值實(shí)驗(yàn)

為了驗(yàn)證模糊C均值聚類小數(shù)據(jù)量法的有效性，以典型的Lorenz系統(tǒng)[11]為例，選兩組不同的系數(shù)進(jìn)行數(shù)值實(shí)驗(yàn).Lorenz系統(tǒng)為

(19)

取4個(gè)不同的初始點(diǎn)A1=[0.046 2，0.097 1，0.823 5]，A2=[0.694 8，0.317 1，0.950 2]，A3=[1，4，10]，A4=[0.119 0，0.498 4，0.959 7]，取a=16，b=4，c=45.92，此時(shí)LLE理論值為1.501 5.取積分步長(zhǎng)為0.01，演化點(diǎn)3 000個(gè)，用小數(shù)據(jù)量法和模糊C均值聚類小數(shù)據(jù)量法(重構(gòu)維數(shù)為5，時(shí)間延遲為10，平均周期為40)計(jì)算LLE所得結(jié)果與誤差率見表1，小數(shù)據(jù)量法簡(jiǎn)稱AL1，模糊C均值小數(shù)據(jù)量法簡(jiǎn)稱AL2.

取4個(gè)不同的初始點(diǎn)A1=[0.046 2，0.097 1，0.823 5]，A2=[0.694 8，0.317 1，0.950 2]，A3=[1，4，10]，A4=[0.119 0，0.498 4，0.959 7]，取a=16，b=4，c=45.92，此時(shí)LLE理論值為1.501 5.取積分步長(zhǎng)為0.01，演化點(diǎn)5 000個(gè)，用小數(shù)據(jù)量法和模糊C均值聚類小數(shù)據(jù)量法(重構(gòu)維數(shù)為5，時(shí)間延遲為10，平均周期為40)計(jì)算LLE所得結(jié)果與誤差率見表2.

取初始點(diǎn)A1=[0.046 2，0.097 1，0.823 5]，取a=16，b=4，c=45.92，此時(shí)LLE理論值為1.501 5.取積分步長(zhǎng)為0.01，分別演化3 000個(gè)點(diǎn)與5 000個(gè)點(diǎn)，用模糊C均值聚類小數(shù)據(jù)量法(重構(gòu)維數(shù)為5，時(shí)間延遲為10，平均周期為40)每次分別演化i，2i，3i，4i(i的取值見上文)個(gè)離散時(shí)間步長(zhǎng)，LLE計(jì)算結(jié)果與運(yùn)行時(shí)間見表3.

取初始點(diǎn)A1=[0.046 2，0.097 1，0.823 5]，取a=10，b=8/3，c=28，此時(shí)LLE理論值為0.881 0.取積分步長(zhǎng)為0.01，分別演化3 000個(gè)點(diǎn)與5 000個(gè)點(diǎn)，用模糊C均值聚類小數(shù)據(jù)量法(重構(gòu)維數(shù)為5，時(shí)間延遲為25，平均周期為40)每次分別演化i，2i，3i，4i(i的取值見上文)個(gè)離散時(shí)間步長(zhǎng)，LLE計(jì)算結(jié)果與運(yùn)行時(shí)間見表4.

從表1與表2結(jié)果可以看出，對(duì)比于小數(shù)據(jù)量法，文中提出的模糊C均值小數(shù)據(jù)量法計(jì)算LLE可以提高計(jì)算精度；從表3與表4結(jié)果可以看出，模糊C均值小數(shù)據(jù)量法計(jì)算LLE時(shí)隨著每次迭代的離散時(shí)間步長(zhǎng)的增加，計(jì)算時(shí)間將縮短，但每次迭代的離散時(shí)間步長(zhǎng)不宜過大.實(shí)驗(yàn)結(jié)果驗(yàn)證了算法是有效可行的.

表1 3 000個(gè)演化點(diǎn)時(shí)小數(shù)據(jù)量法與模糊C均值小數(shù)據(jù)量算法的計(jì)算結(jié)果對(duì)比

表2 5 000個(gè)演化點(diǎn)時(shí)小數(shù)據(jù)量法與模糊C均值小數(shù)據(jù)量法的計(jì)算結(jié)果對(duì)比

表3 模糊C均值小數(shù)據(jù)量法改變步長(zhǎng)的計(jì)算結(jié)果1(LLE理論值為1.501 5)

表4 模糊C均值小數(shù)據(jù)量法改變步長(zhǎng)的計(jì)算結(jié)果2(LLE理論值為0.881 0)

4 結(jié)論

將模糊C均值聚類算法用于小數(shù)據(jù)量法線性區(qū)域的選擇而提出的模糊C均值小數(shù)據(jù)量法來計(jì)算LLE，通過選出小數(shù)據(jù)量法的線性區(qū)域，可以大大提高小數(shù)據(jù)量法的計(jì)算精度.實(shí)驗(yàn)表明，當(dāng)每次迭代的離散時(shí)間步長(zhǎng)增加時(shí)，計(jì)算時(shí)間將減少.接下來可以研究離散步長(zhǎng)的選取對(duì)于實(shí)驗(yàn)精度和效率的影響，從而得到更為完善的模糊C均值小數(shù)據(jù)量法.