9512.net

# 不光滑混沌系统的Lyapunov指数计算

Computation of Lyapunov Exponents of Non-smooth Chaotic Systems

by Wang Xinxin

Supervisor: Associate Professor Zhang Xilin

Tsinghua University June 2006

1. 对混沌理论和 Lyapunov 指数相关信息进行研究； 2. 仔细研究相空间重构理论和求 Lyapunov 指数的算法； 3. 确定算法，并进行编程； 4. 求解 Chua 电路和振动系统的 Lyapunov 指数，并分析是否混沌； 5. 总结全文。

I

II

Abstract

Computation of Lyapunov Exponents of Non-smooth Chaotic Systems

Abstract
Chaos theory is considered as the one of three scientific revolutions in 20th century. From the time that chaotic theory appeared to now, as the rapid development of computer technology, more and more scholars concerned about the chaos theory, and chaos theory has been so tremendous development. Where, a very important indicator to judge whether the system is chaotic is the Lyapunov exponent. It can be easier to judge whether the system is chaotic. First, this dissertation introduces chaos theory and the Lyapunov exponents, including their definition, history, and diagnosis method. Then, this article introduces the reconstruction of phase space and Lyapunov exponents detailedly, including their definition, character and algorithms. Algorithms used in this dissertation are Wolf and CC method (CC method is used to calculate embedding dimension and delay time, Wolf is used to calculate the maximum Lyapunov exponent), and these algorithms are programmed. And then, this dissertation calculates Lyapunov exponents of Chua circuit and vibration system. Then changing the parameters, calculate the Lyapunov exponents, judge that systems with different parameters is chaotic. Illustrates the sensitivity of that chaotic systems dependent on the parameters and initial value. Finally, conclude this dissertation, prospect this thesis.

Keywords: Lyapunov exponent, phase space reconstruction, Chua circuit, chaos, Wolf method

III

3.2 关于 Lyapunov 指数的性质 ...................................................................................... - 21 3.3 计算 Lyapunov 指数的方法 ...................................................................................... - 23 3.4 本章小结 .................................................................................................................... - 24 第四章 具体算法及实现 ................................................................................................... - 25 4.1 微分方程求解 ............................................................................................................ - 25 4.2 C-C 算法的具体实现 ................................................................................................ - 26 4.2.1 子程序设计： .................................................................................................... - 27 4.2.2 主程序设计 ........................................................................................................ - 29 4.3 Wolf 算法的具体实现 ............................................................................................... - 31 4.4 本章小结 .................................................................................................................... - 32 第五章不光滑系统的 Lyapunov 指数和分析.................................................................. - 33 指数和分析 5.1 Chua 电路及其 Lyapunov 指数计算 ........................................................................ - 33 5.1.1 Chua 电路介绍................................................................................................... - 33 5.2.2 Chua 电路的化简............................................................................................... - 36 5.3.3 不同参数下 Chua 电路的混沌性...................................................................... - 36 5.2 振动系统的 Lyapunov 指数计算 .............................................................................. - 40 5.2.1 振动系统的介绍和化简 .................................................................................... - 40 5.2.2 振动系统的 Lyapunov 指数计算 ...................................................................... - 42 5.3 本章小结 .................................................................................................................... - 44 第六章 总结与展望 ........................................................................................................... - 45 参考文献 ............................................................................................................................. - 47 致谢 ..................................................................................................................................... - 51 -

V

1.1 背景和意义

1.2 混沌理论和相关知识
1.2.1 混沌的定义

d ≤ a < b < c 或 d ≥ a > b > c ， f 映射按 Li.Yorke 定义。

-1-

Sf ( x) =

f m ( x) 3 f n ( x) 2 ? ( ) <0 f ′( x) 2 f ′( x)

(1.1)

(1) 敏感地依赖于初值：存在 δ > 0 ，对于任意的 ε > 0 和任意 x∈V，在 x 的 ε 邻域

(2) 拓扑传递性：对于 V 上任意～对开集 X，Y，存在 k>0，使 fk ( x) ∩ Y ≠ φ 。 (3) f 的周期点集在 V 中稠密。

1.2.2 混沌发展史

[12]

， 也有学者研究大型输变电电力系统混沌模型[13]， 也有文章报道应用混沌动力学研究

1.2.3 混沌的基本特征

-6-

1.2.4 混沌的判别方法

1.3 Lyapunov 指数
1.3.1 Lyapunov 指数的现状

1.3.2 Lyapunov 指数的性质
Lyapunov 指数是衡量系统动力学特性的一个重要定量指标,它表征了系统在相空间 中相邻轨道间收敛或发散的平均指数率。对于系统是否存在动力学混沌,可以从最大 Lyapunov 指数是否大于零非常直观的判断出来：一个正的 Lyapunov 指数，意味着在系 统相空间中，无论初始两条轨线的间距多么小，其差别都会随着时间的演化而成指数率 的增加以致达到无法预测,这就是混沌现象。 Lyapunov 指数的和可以表示椭球体积的变化率。 对于汉米尔顿系统， 它的 Lyapunov 指数的和等于零；对耗散系统,Lyapunov 指数的和为负。若耗散系统的吸引子是一个不 动点，那么它的 Lyapunov 指数一般都是小于零的。如果是一个简单的 m 维流形（m=1 或 m=2 分别为一个曲线或一个面）那么， m 个 Lyapunov 指数等于零,其余的 Lyapunov ， 前 指数小于零。总之，不管所求系统是不是耗散的，只要 λ1>0 就会出现混沌。 微分动力系统 Lyapunov 指数的性质： 对于一维(单变量)情形，吸引子只可能是不动点(稳定定态)，此时 λ 是负的。对于 二维情形，吸引子或者是不动点或者是极限环。对于不动点，任意方向的 δxi，都要收 缩，故这时两个 Lyapunov 指数都应该是负的，即对于不动点，(λ1,λ2)=(-,-)。至于极限环， -8-

1.4 本文内容介绍

P-范数，Wolf 方法等等，其中对 Wolf 方法进行了详细说明。 第四章是对于本文中涉及到的编程和绘图方法进行介绍。 第五章是几个不光滑混沌系统的 Lyapunov 指数计算实例， 其中包括著名的 Chua 电 路和振动系统。通过实例深入研究 Lyapunov 指数与混沌系统的关系。 第六章是对于本课题的总结和展望。

- 10 -

2.1 概述

2.2 相空间重构相关理论
2.2.1 相空间重构理论的概述

2.2.2 Taknes 定理

Takens 的嵌入定理是在非常一般的情况下讨论的。定理中当 ? 是某些特定函数， M 是非紧流形时，定理同样成立。函数 g 实质上就是 D 维流形 L 中人们所能观测到的

g ( x(t )) = g (t ) , g (? ( x(t ))) = g ( x(t ? τ )) = g (t ? τ ),? , g (? 2 D ( x(t ))) = g (t ? 2 Dτ ) 。这样嵌入定

2.3 重构相空间参数的求取

2.3.1 求取时间延迟 τ 的方法

1、互信息(mutualinformation)法[41]； 2、其它信息论的方法，如冗余度(redundancy)法、信息熵等方法； 3、自相关法和复自相关法[41]； 4、预报效果法； 5、真实矢量场法； 6、波动积(wavering product)法； 7、填充因子(fill factor)法； 8、累积局部变形(integral local deformation)法；

- 13 -

9、简单轨道扩张(simple spreading of trajectories)法； 10、奇异值分数(SVF)法和 Jacobi 矩阵行列式法。

2.3.2 求取嵌入维数 m 的方法

2.3.3 同时求取嵌入维数 m 和时间延迟 τ 的方法（C-C 法）

Kugiumtizs 提出时间延迟 τ d 的选取不应该独立于嵌入维数 m，而应该依赖于延迟时间窗

τ s 值时间序列的采样间隔， τ d = tτ s 指时间序列的延迟， τ w = (m ? 1)τ d 指延迟时间窗

X i = ( xi , xi +τ ,? , xi + ( m ?1)τ ), X i + ∈ R m

C (m, N , r , t ) = 2 ∑ θ (r ? dij ), r > 0 M ( M ? 1) 1≤i ≤ j ≤ M

(2.1)

(2.2)

D(m, t ) = lim

(2.3)

log C (m, r , t ) r →∞ log r

(2.4)

C (m, r , t ) = lim C (m, N , r , t )
N →∞

(2.5)

D(m, t ) = log C (m, N , r , t ) log r (2.6)

{x1 , xt +1 , x2t +1 ,?} {x2 , xt + 2 , x2t + 2 ,?}

??

(2.7)

{xt , x2 t , x3t ,?}

N = tl , l =

N 是长度。对于一般的时间序列，将其分成 t 个不相交的子序列，然后定义每 t

S (m, N , r , t ) = 1 t ? N N ? ∑ ?Cs (m, t , r, t ) ? Csm (1, t , r , t ) ? t s =1 ? ? (2.8)

S ( m, r , t ) =

1 t ∑ ?Cs (m, r , t ) ? Csm (1, r , t ) ? , (m = 2,3,?) ? t s =1 ?

(2.9)

S (m, r , t ) 一般不等于零。这样，局部最大时间间隔可以取 S (m, r , t ) 的零点或对所有的半

?S (m, t ) = max{S (m, rj , t )} ? min{S (m, rj , t )}

(2.10)

?S (m, t ) 的最小值。但是， S (m, r , t ) 的零点对所有 m, r 应几乎相等； ?S (m, t ) 的最小值对

σ
2

≤ r ≤ 2σ , N ≥ 500 ，渐近分布可以通过有限序列

iσ , i = 1, 2,3, 4 。计算 2

- 16 -

S (t ) =

1 5 4 ∑∑ S (m, rj , t ) 16 m = 2 j =1
1 5 ∑ ?S (m, t ) 4 m=2

(2.11)

? S (t ) =

(2.12) (2.13)

Scor (t ) = ? S (t ) + S (t )

(1) 计算给定时间序列的标准差 σ ，选取 N=3000。 (2) 编程在计算机上计算下列三个量

S (t ) =

1 5 4 ∑∑ S (m, rj , t ) 16 m = 2 j =1
1 5 ∑ ?S (m, t ) 4 m= 2

(2.14)

? S (t ) =

(2.15) (2.16)

Scor (t ) = ? S (t ) + S (t )

S (m, rj , t ) = 1 t ∑ ?Cs (m, rj , t ) ? Csm (1, rj , t ) ? , (m = 2, 3, 4,5) ? t s =1 ? (2.17) (2.18)

?S (m, t ) = max{S (m, rj , t )} ? min{S (m, rj , t )} (3) 依据上述计算结果画图 (a) ?S (m, t ) → t (0 ≤ t ≤ 200) ， ? S (t ) 的第一个极小值 t，对应延迟 τ d = tτ s 。 (b) ? S (t ) → t (0 ≤ t ≤ 200) ， ? S (t ) 的第一个极小值 t 对应延迟 τ d = tτ s 。 S (t ) → t (0 ≤ t ≤ 200) ， S (t ) 的第一个零点 t 对应时间延迟 τ d = tτ s 。 (c) Scor (t ) → t (0 ≤ t ≤ 200) ，最小值 t 对应时间窗口 τ w = tτ s 。
v. s .
v.s. v. s . v. s .

C-C 方法的特点是： 1、容易操作；2、计算量小；3、对小数据组可靠；4、效果和互信息方法一致；5、

。 具有较强的抗燥声能力（30%以下） - 17 -

2.4 本章小结

- 18 -

3.1 Lyapunov 指数的定义

(3.1)

dF > 1 ，映射使两点分开 dx dF < 1 ，映射使两点靠拢 dx

- 19 -

?F2

?F1

?x1

?x2

x

ε eλ
x0

x0 + ε

F ( x0 )

F ( x0 + ε )

εe

F 2 ( x0 )

F 2 ( x0 + ε )

ε e nλ
n次迭代
F n ( x0 )

F n ( x0 + ε )

ε e nλ ( x ) = F n ( x0 + ε ) ? F n ( x0 )
0

(3.2)

λ = lim ∑ ln
n →∞ i =0

n

dF dx

(3.3)
x = xi

λ 称为 Lyapunov 指数(Lyapunov exponent)，它表示大量次数迭代中，平均每次迭代所引

λi = lim lim

1 ε i (t ) ln t →∞ ε (0) → 0 t ε (0)

(3.4)

ε (0)

ε (0)

ε1

ε3

ε2

3.2 关于 Lyapunov 指数的性质

δ xi 都要收缩， 故这时两个 Lyapunov 指数都是负的：λ1 < 0 ，λ2 < 0 。 至于极限环， δ xi 取

λ1 + λ2 + ? + λk 决定。而 ∑ λi 则决定整个 n 维超球体的体积收缩的快慢。由此可见，对
i =1

n

∑ λ ?= 0，对于保守系统
i =1 i

n

?< 0，对所有耗散系统 ?

(2) 稳定定态和周期运动(以及准周期运动)不可能有正的 Lyapunov 指数。稳定定态 的 λi 都是负的；周期运动的 λ1 = 0 ，其余 λi 也都是负的； (3) 混沌运动至少有一个正的 Lyapunov 指数。反之，如果计算得知系统至少有一个 正的 Lyapunov 指数。 ，则可肯定系统作混沌运动。 现在可以用最大李雅普诺夫指数的取值定量描述混沌现象了， 就是用最大 Lyapunov 指数 λ1 的正负来区别，而正的 λ1 的大小则可看作是混沌强弱程度的定量标志。

3.3 计算 Lyapunov 指数的方法

λ=

1 tM ? tOi

∑ ln L
i =0

M

Li′
i

(3.5)

λ1 + λ2 =

1 tM ? tOi

∑ ln A
i =0

M

Ai′
i

(3.6)

L′(t1 )

L′(t3 ) L′(t2 )

θ1
L(t0 )
t0 t1

L(t1 )

L(t2 )
t2

θ2

t3

M ′(t2 )

M ′(t1 )

M ′(t3 )

t3

M (t1 )
M (t0 ) t0
t1

M (t2 )
t2

3.4 本章小结

- 24 -

4.1 微分方程求解

? y ( n ) = f (t , y, y,? , y ( n ?1) )

(4.1)

? 其中输出变量 y (t ) 的各阶导数初始值为 y (0), y (0),? , y ( n ?1) (0) ，则可以选择一组状态变量
? 来进行替换 x1 = y, x2 = y,? , xn = y ( n ?1) ，于是原方程就变换为：

? x1 = x2 ? ? ? x2 = x3 ? ? ? ? ? xn = f (t , x1 , x2 ,? , xn ) ?

(4.2)

? 上式的初值为 x1 (0) = y (0) ， x2 (0) = y (0) ，…， xn (0) = y ( n ?1) (0) 。于是，变换后就是比较

? ? x1 (t ) = ? β x1 (t ) + x2 (t ) x3 (t ) ? ? ? x2 (t ) = ? ρ x2 (t ) + ρ x3 (t ) ? x (t ) = ? x (t ) x (t ) + σ x (t ) ? x (t ) ? ?3 1 2 2 3

(4.3)

4.2 C-C 算法的具体实现

C-C主程序

Heaviside.m

Reconstitution.m

Disjoint.m

Guanlianjifen.m

4.2.1 子程序设计：
C-C 算法中要用到四个子函数： (1) Heaviside.m：用来求解 Heaviside 函数的值； (2) Reconstitution.m：用来进行相空间重构； (3) Disjoint.m：用来将时间序列分拆为 t 个不相关的时间序列； (4) Guanlianjifen.m：用来计算时间序列的关联积分。 程序说明： 1、Heaviside.m：

θ ( x) = ?

?0, x < 0 ?1, x ≥ 0

(4.4)

%该函数用来重构相空间 % data为输入时间序列 % N为时间序列长度 % m为嵌入空间维数 % t为时间延迟 % X为输出，是m*n维矩阵 M=N-(m-1)*t;%相空间中点的个数 for j=1:M for i=1:m X(i,j)=data((i-1)*t+j); end end 3、disjoint.m：将时间序列分成 t 个不相交的时间序列，长度 l = N / t
{x1 , xt +1 , x2t +1 ,?} {x2 , xt + 2 , x2t + 2 ,?}

%相空间重构

???
{xt , x2t , x3t ,?}

function data_d=disjoint(data,N,t) %这个程序是用来将时间序列分解成不相关的时间序列 %data:输入时间序列 %N:时间序列的长度 %t:延迟时间 %data_d:返回的不相关时间序列 for i=1:t for j=1:(N/t) data_d(i,j)=data(i+(j-1)*t); end end 4、Guanlianjifen.m：

2 ∑θ (r ? d ij ), r > 0 M ( M ? 1) 1≤i≤ j ≤ M

(4.5)

- 28 -

θ ( x) = ?

?0, x < 0 ?1, x ≥ 0

(4.6)

(1) 输入数据 X：X 为根据给定的时间延迟 t 和嵌入维数 m 进行相空间重构后得到的
m × M 的矩阵；

(2) 对相空间中所有点求距离，并保存在数组 d(i,j)中； (3) 调用 Heaviside 函数计算所有的 Heaviside 函数值并求和得到 H_sum (4) 计算出关联积分的值，并输出 Jifen_out。 function Jifen_out=guanlianjifen(X,M,r) %这个程序是用来计算关联积分 %Jifen_out:关联积分的值 %X:重构的相空间,M是一个m*M矩阵 %m:嵌入维数 %M:M是在m维空间中嵌入点的个数 %r:Heaviside函数的半径,范围sigma/2<r<2sigma %H_sum：计算所有Heaviside返回值的和 H_sum=0; for i=1:M % fprintf('%d/%d\n',i,M); for j=i+1:M d=norm((X(:,i)-X(:,j)),inf);%计算两向量之间距离的范数 sita=heaviside(r,d);%计算heaviside函数的值 H_sum=H_sum+sita;%对sita求和 end end Jifen_out=2*H_sum/(M*(M-1));%关联积分的值

4.2.2 主程序设计
C-C 中的主程序程序步骤： (1) 读入数据； - 29 -

(2) 计算时间序列的标准差； (3) 让 r 从 sigma/2 到 2sigma 变化，m 从 2 到 5 变化，t 从 1 到 200 变化； (4) 调用 disjoint 函数将时间序列分拆成 t 个不相交的子列，并调用 guanlianjifen 函数计 算 C(1,N,r,t); (5) 调用 reconstitution 函数对子列进行相空间重构，并调用 guanlianjifen 函数计算 C(m,N,r,t); (6) 计算 C(m,N,r,t)- C(1,N,r,t)m，并对 t 求和得到 s_t3； (7) 根据公式 S (m, r j , t ) =
s_t1； (8) 将 s_t1 赋给 s_t0(m)，并对 m 求和，得到 s_t； (9) 根据公式 S (t ) = 1 t ∑ [Cs (m, rj , t ) ? Cxm (1, rj , t )] 计算得到 s_t2(j)，并对 r j 求和得到 t s=1

1 5 4 ∑∑ S (m, r j , t ) 求得 s(t)； 16 m= 2 j =1

(10) 同 时 利 用 已 求 出 的 s(t) ， 根 据 公 式 ?S (m, t ) = max{S (m, r j , t )} ? min{S (m, r j , t )} 和

S cor (t ) = ? S (t ) + S (t ) 分别求出相应的 delt_s(t)和 s_cor(t);
(11) 根据求得的结果作图。

figure; subplot(311); plot(1:maxLags,S_mean);grid;title('S mean'); subplot(312); plot(1:maxLags,delta_S_mean);grid;title('delta S mean'); subplot(313); plot(1:maxLags,S_cor);grid;title('S cor');

- 30 -

4.3 Wolf 算法的具体实现
Wolf 方法提取 Lyapunov 指数： (1) 应用时间序列选取时间延迟 τ ，根据观测数据样本总数 N，构造 m 维相空间的 新序列，相点数为 M： M = N ? (m ? 1)τ ；
(2) 以初始相点 X 0 为基点，在点集 { X i } 的其余相点中选取与 X 0 最近的点 X j 为端

(3) 取时间步长为 k，t1 = t0 + k ，初始向量沿轨线向前演化得到一新向量，其相应基

λ1 = log 2

1 k

L(t1 ) ； L(t0 )

(4) 如此继续直至所有相点， 然后取各指数增长率的平均值为最大 Lyapunov 指数估

1 M

∑ k log
i =1

M

1

2

L(ti ) ； L(ti ?1 )

(5) 依次增加嵌入维数 m，重复(2)、(3)、(4)过程直至 Lyapunov 指数估计值随 m 变

Lyapunov 指数的计算结果影响不大，加之现在没有有效的计算 P 的方法，所以 P 看图

- 31 -

4.4 本章小结

- 32 -

5.1 Chua 电路及其 Lyapunov 指数计算
5.1.1 Chua 电路介绍
20 世纪 80 年代初，随着非线性动力学系统研究的蓬勃发展，电路的混沌现象越来 越被更多的学者关注，并设计了多种电路来模拟，而其中最为经典的就是 Chua 电路。 在非线性理论中，一个三维自治系统能够产生混沌至少有一个稳定的不动点和两个 不稳定的不动点，具体到电路而言，该电路必须满足三个条件： (1) 至少一个非线性元件； (2) 至少一个有源阻抗； (3) 至少三个储能元件。 而 Chua 电路是符合上述条件的最简单的电路，如图 5.1 所示。

R

L

C2

C1

Nr

Chua 是 1983 年 Leon. O. Chua 提出的并且通过数值计算证明可以产生混沌的第一个 物理系统，通过实际电子电路实验，也被观测到的混沌现象所确认。Chua 电路包括一 个线性电感 L、两个线性电容 Cl、C2、一个线性电阻 R 和一个负阻器件 Nr。Chua 电路 的动力学行为可以用一个三阶自治方程——方程组(5.1)，来描述。 - 33 -

dVc1 Vc 2 ? Vc1 1 = ? f (Vc1 ) dt RC1 C1 dVc 2 Vc1 ? Vc 2 = +i dt RC2 di 1 = ? Vc 2 dt L 上面式子中 Vc1 ，Vc 2 ，i 分别为电容 C2，C2 两端的电压和流过电感 L 的电流, f (Vc1 ) 是描述负阻器件 Nr 的 i-Vc1 特性曲线函数。
Nr 一般被称作 Chua 二极管，是一个非线性元件，可以用运算放大器、二极管或三

(5.1)

VDD 15V

R3 22k
8 3 1 2 4

R6 220 U1A
3 1 2 8

U2A

TL082CD R5 220 R4 2.2k

4

TL082CD

R2 22k R1 3.3k

VSS -15V

?VR V0 ? VR ? = R2 ? R1 ?V0 ? VR = ?iR R3 ?

R2 R2 VR ，则其等效阻抗为： G1 = ? R1 ? R3 R1 ? R3

- 34 -

(2) 当 V0 输出达到饱和值 E 时， V0 为常值，此时 iR = 抗： G 2 =

1 E VR ? ，则得到其等效阻 R3 R3

1 ，E 为运算放大器的输出饱和电压，跟运算放大器的工作电压有关。 R3 R2 1 E VR ， iR = VR ? 两式同时成 R1 ? R3 R3 R3

(3) 在预算放大器饱和输出的临界点， iR = ?

V1 =

1 E R2 VR ? =? VR ， 得 到 临 界 点 的 电 压 即 图 5.3 中 的 转 折 点 电 压 R3 R3 R1 ? R3

1 E 。同理，与其并联的另一个电路的各参数求解式子如下： R2 1+ R1 G1′ = ? R5 1 1 ， G 2′ = ，V 2 = E R5 R 4 ? R6 R6 1+ R4

Ga = G1 + G1′ = ? R2 R5 ? R1 ? R3 R 4 ? R6 1 R5 ? R3 R 4 ? R 6 1 1 + R3 R 6

Gb = G 2 + G1′ =

Gc = G 2 + G 2′ =

?GcVC1 + (Gc ? Gb)V2 + (Gb ? Ga )V1 VC1 < ?V2 ?GbV + (Gb ? Ga )V ?V2 ≤ VC1 < ?V1 1 C1 ? ? f (VC1 ) = ?GaVC1 VC1 ≤ V1 ?GbV + (Ga ? Gb)V V1 ≤ VC1 < V2 1 C1 ? ?GcVC1 + (Gc ? Gb)V2 + (Gb ? Ga )V1 VC1 > V2 ?

(5.2)

1 f (V 1) = GbV1 + (Ga ? Gb) × [ V1 + E ? V1 ? E ] 2

(5.3)

- 35 -

XSC1
Ext T rig + _ A + _ + B _

VDD 12V

R 1500 L 18.68mH C2 100nF C1 10nF

R3 39k
8 3 1 2 4

R6 220 U1A
3 1 2 8

U2A

TL082CD R5 220 R4 2.2k

4

TL082CD

R2 39k R1 3.3k

VSS -12V

5.2.2Chua 电路的化简

? dx ? d τ = α ( y ? x ? f ( x )) ? ? dy = x? y+ z ? ? dτ ? dz ? dτ = ? β y ?

(5.4)

(5.5)

1 式(5.5)中， f ( x) 为 f ( x) = bx + (a ? b)( x + 1 ? x ? 1) ，它是一个三段的折线方程，a 和 b 2

5.3.3 不同参数下 Chua 电路的混沌性

R=2100 时， a=-1.5909,b=-0.9007, α =10, β =23.6081。 嵌入维数 m=2， 延迟时间 τ=80。 计算 Lyapunov 指数为-0.0048 - 37 -

R=2000 时， a=-1.5152,b=-0.8578, α =10, β =21.4133 嵌入维数 m=3， 延迟时间 τ=24。 计算 Lyapunov 指数为-0.0040

R=1900 时，a=-1.4394,b=-0.8149, α =10, β =19.3255

(a) R=1800 的相图

(b) R=1700 的相图

(c) R=1600 的相图

R=1800

- 39 -

R 由 1594

5.2 振动系统的 Lyapunov 指数计算
5.2.1 振动系统的介绍和化简

a0 k1 m f
1

a k2 m0 f
2

?? ? mx + F ( x, x) + Q ( x) = m0ω 2 r sin(ωt )
? 式(5.6)中， F ( x, x) , Q ( x) 分别表示振动系统的阻尼力和弹性恢复力，其中：

(5.6)

? ?( f + f ) x ? F ( x, x ) = ? 1 2 ? f1 x ?(k1 + k2 ) x Q ( x) = ? ?k1 x + k1 (a0 ? a)

x > ?a x ≤ ?a x > ?a x ≤ ?a
(5.7) (5.8)

m ——振动质量体质量
k1 ——软弹簧刚度 f1 ——软弹簧阻尼

m0 ——偏心块质量 k2 ——硬弹簧刚度 f 2 ——硬弹簧阻尼
r ——偏心量
? x ——振动质量体速度 a0 ? a ——静平衡下软弹簧压缩量

ω ——激振角频率
x ——振动质量体位移
?? ——振动质量体加速度 x

a ——静平衡下硬弹簧压缩量

- 41 -

-a 速度 -a f1

5.2.2 振动系统的 Lyapunov 指数计算

? 令 y (1) = x, y (2) = x ，则：

? my (1) + F ( y (1), y (2)) + Q ( y (2)) = m0ω 2 r sin(ωt )

(5.9) (5.10)

?( f + f ) y (1) F ( y (1), y (2)) = ? 1 2 ? f1 y (2)

y (2) > a y (2) ≤ ? a

y (2) > ? a ?(k1 + k2 ) y (2) Q( y (2)) = ? ?k1 y (2) + k1 (a0 ? a ) y (2) ≤ ? a

(5.11)

? dy (2) ? dt = y (1) ? ? 2 ? dy (1) = m0ω r sin(ωt ) ? ( f1 + f 2 ) y (1) ? (k1 + k2 ) y (2) ? dt m ? y(1)≤-a 时， ? dy (2) ? dt = y (1) ? ? 2 ? dy (1) = m0ω r sin(ωt ) ? f1 y (2) ? (k1 y (2) + k1 (a0 ? a )) ? dt m ?

(5.12)

(5.13)

m=400，m0=25，k1=43700，k2=900000，f1=500，f2=50，ω =68.56，r=0.126，a0=5，a=10。

- 42 -

(m0*(w.^2)*r*sin(w*t)-(f1+f2)*x(2)+(k1+k2)*x(1))/m]; else xdot=[x(2); (m0*(w.^2)*r*sin(w*t)-f1*x(1)-k1*x(1)-k1*(a0-a))/m]; end 有了上面程序，对振动系统求解，然后作图。得到两列时间序列，并且得到系统的 时域图和相图。然后由 C-C 法求得系统的 S、?S 和 Scor 图，然后分析得到了系统的嵌入 维数 m 为 4， 延迟时间 τ 为 42。 然后再计算 Lyapunov 指数， 得到的 0.0280。 最大 Lyapunov 指数为正，证明系统在上述参数下是混沌系统。 %%对方程(5.13)编程

- 43 -

5.3 本章小结

- 44 -

- 45 -

- 46 -

1. Lorenz E N. Deterministic non-periodic flow[J]. J. Atmos. Sci., 1963, 20: 130-141. 2. Henon M. A two-dimensional mapping with a strange attractor[J]. Comm. Math. Phys, 1976, 50(1): 69-77. 3. 张济忠. 分形[M]. 北京：清华大学出版社，1995． 4. 匡锦瑜，邓昆，黄荣怀. 利用时空混沌同步进行数字加密通信[J]. 物理学报，2001， 50(10)：1856-1861. 5. 谢鲲， 雷敏， 冯正进. 一种超混沌系统的加密特性分析物理学报[J]. 物理学报， 2005， 54(3)：1265-1272. 6. 方锦清，王中生，陈关荣. 若干束混沌的非线性控制方法的稳定性分析[J]. 中国原子 能科学研究院年报，2004，2004：106-107. 7. 方锦清. 高科技发展的新领域——驾驭混沌的应用前景[J]. 科学导报， 2003， 3-6． 10： 8. 王林翔，陈鹰，路甬祥. 二次流实现流体混沌混合的数值研究[J]. 化工学报，1999， 50(4)：449-455. 9. 郑会永， 肖田元， 韩向利等. 心电信号的混沌分形特性分析[J]. 清华大学学报， 1999， 39(9)：103-105. 10. 刘耀宗. 碰摩转子混沌振动识别与控制技术研究[D]. 长沙：国防科技大学，2001. 11. 何斌，杨灿军，陈鹰. 混沌周期解提高测量灵敏度算法及抗干扰分析[J]. 电子学报， 2003，3l(1)：68-70. 12. 蒋海峰， 王鼎媛， 张伸义. 短时交通流的非线性动力学特性[J]. 中国公路学报， 2008， 21(3)：91-96. 13. 王东风，韩璞，于朝辉. 电力系统中的混沌研究与混沌应用[J]. 电力科学与工程， 2003，3：74-78. 14. 郭建青，李彦，王洪胜等. 利用混沌优化算法确定河流水质模型参数[J]. 水力发电学 报，2004，23(4)：92-96. 15. 修日晨， 张自历， 刘爱菊. 海洋激流的观测实验及分析讨论[J]. 海洋学报， 2004， 26(2)： 118-124. 16. 屈元举，谢芳森，沈艳菲. 一个新混沌系统的动力学特性研究和电路仿真[J]. 江西科 学，2010，28(1)：1-4. - 47 -

17. 李银山，李欣业，刘波. 分岔混沌非线性振动及其在工程中的应用[J]. 河北工业大学 学报，2004，33(2)：96-103. 18. 张淑芹， 姚翠珍. Euler-Bernoulli 梁在自激励边界反馈下的动力行为[J]. 数学的实践与 认识，2004，34(7)：67-71. 19. 胡春林，彭华中，程昌钧. 材料及结构参数对基桩非线性纵向振动的影响[J]. 武汉理 工大学学报，2006，28(4)：71-74. 20. 金基铎，邹光胜，闻邦椿. 参与强激励联合作用下输流管的分岔和混沌行为研究[J]. 应用力学学报，2005，22(1)：111-113. 21. 欧阳茹荃，朱继懋. 船舶横摇运动的非线性振动与混沌[J]. 水动力学研究与进展， 1999，14(3)：334-340. 22. 姜荣俊，朱石坚，何琳. 混沌隔振装置的设计[J]. 海军工程大学学报，2003，15(1)： 35-37. 23. 韩正铜，张永忠，黄民. 磨削过程振动特征的实验研究[J]. 制造技术与机床，2004， 2：63-65. 24. 王跃科， 林嘉宇， 黄芝平等. 语音信号非线性分析与处理[J]. 通信技术， 2010， 108(1)： 61-65. 25. 李亚安，贾雪松，孙进才. 基于局部投影理论的水声信号降噪处理研究[J]. 西北工业 大学学报，2005，23(2)：147-151. 26. Frey M, Simi E. Noise-induced chaos and phase space flux [J]. Physica D, 1993, 63: 321-340. 27. 宋春云，丁士圻，雷亚辉. 基于相空间重构理论的水中混响特性提取[J]. 哈尔滨工程 大学学报，2007，28(5)：515-518. 28. 丁旺才，谢建华. 碰撞振动系统分岔与混沌研究进展[J]. 力学进展，2005，35(4)： 513-524. 29. 章新华， 张晓明， 林良骥. 船舶辐射噪声的混沌现象研究[J]. 声学学报， 1998.3， 23(2)： 134-140. 30. 杨杰，卢凌. 海杂波的混沌分析[J]. 交通与计算机，1998，16(4)：7-10. 31. 陈亮，张雄伟. 语音信号非线性特征的研究[J]. 解放军理工大学学报，2000，l(2)： 11-17. 32. 王兴元，朱伟勇. 利用心电的李氏数探讨混沌在生物进化中的作用[J]. 清华大学学 报，1999，20(3)：247-249. - 48 -

33. 谢文会，陈予恕. 利用 Lyapunov 指数分析一种非稳态油膜轴承单盘柔性转子的分岔 和混沌特性[J]. 非线性动力学学报，2002，9(3)：145-147. 34. 杨文平，陈国定，石博强等. 基于李雅普诺夫指数的汽车发动机故障诊断研究[J]. 振 动工程学报，2002，15(4)：476-478. 35. 侯荣涛，褚祥忐，仟立义. 基于 Lyapunov 指数相轨迹的发电机组故障诊断[J]. 清华 大学学报，2004，25(5)：478-481. 36. 张来斌，工朝晖，田立柱. 机械设备故障信弓的李雅普诺夫指数识别[J]. 石油矿场机 械，2004，33(2)：5-7. 37. 田韦国. 基于神经网络的 Lyapunov 指数谱的计算[J]. 系统与工程理论与实践，2001， 8：9-13. 38. 陈淑燕，王炜. 基于 Lyapunov 指数的交通量混沌预测方法[J]. 土木工程学报，2004， 37(9)：97-99. 39. 许传华， 任青文， 房定旺. 基于神经网络的混沌时间序列预测[J]. 水文地质工程地质， 2003，1：30-32. 40. Takens F. Determing strange attractors in turbulence[J]. Lecture notes in Math, 1981, 898(2): 361-381. 41. 张淑清， 贾健， 高敏等. 混沌时间序列重构相空间参数选取研究[J]. 物理学报， 2010， 3(59)：2-5. 42. H. S. Kim, R. Eykholt, J. D. Salas. Nonlinear dynamics, delay times, and embedding windows, Physica D, 1999, D127: 48-60. 43. 刘秉正，彭建华. 非线性动力学[M]. 北京：高等教育出版社，2004. 403.408. 44. 梁志珊，王利敏. 基于 Lyapunov 指数的电力系统短期负荷预测[J]. 中国电机工程学 报，1988，18(5)：368-371. 45. A. Wolf. J. B. Swift, H. L. Swinney J. A. Vastano. Determining Lyapunov exponents from a time series[J], Physica 16D, 1985, 285-317

- 49 -

- 50 -

- 51 -