·138· 工程科学学报，第38卷，第1期 (A>0)收缩形成稳定吸引子.根据Lyapunov指数的 从K()=Q'()0()中可以得到 性质可知道，只要有一个Lyapunov指数为正，就可以 0(t)=Q(t)K(t),Q(0)=L (13) 判断该系统处于混沌振荡状态。 方程(12)和方程(13)就确定了关于矩阵Q()的微分 系统(6)的线性变分方程为 方程.由于Q(d)J()Q(t)是斜对称的，由方程(11) (t)=J(t)Y(t),Y(0)=1 (8) 得到 Y∈R3x3,l,是3×3的单位矩阵，J(t)是系统的Jacobi 矩阵 R(t) R.(t） R0=(Q0J0Q(),R(0)=0. rf1)f1)f(1)1 (14) ayu dy2 ay: 这里R(i)=ln(R.(t)),Lyapunov指数的时间演化式 J(t)= (2)f(2) f(2) ay ay2 ays 为A,(t)=R(t)/L,所以得到Lyapunov指数为 f(3)f(3)f(3) R(t） 入，=im=liml:(0,i=l2,3. (15) L ayay, dy; 0 1 利用QR分解算法和MATLAB软件计算得到如图 0 3所示的Lyapunov指数谱，入，>入2>入·可以看到， -1 -K 1 (9) q=12.55时，Lyapunov指数在形态上发生“突变”，在 -B/Y3 0 1 2√J 章节2.1中分析可知，此时发生擦边分岔. (a 由于Y(0)是非奇异的，所以对于t≥0，Y(t)也是 非奇异的.Y(t)可QR因数分解，记为Y(t)=Q(t) R(t),上三角矩阵R(t)的对角元素都是正的.将该 -1 QR分解式带入方程(8)，有 -3 0(t)R()+Q()R()=J()Q()R(), Q(0)R(0)=I· (10) -S 将式(10)左乘Q',右乘R,有 6 Q()o()+R()R-()=Q()J()2(), -7 Q(0)=L3,R(0)=13 (11) 10 2 由于Q()Q(t)=L3,矩阵K(t)=Q(t)0(t)是斜对 角矩阵，并且由于矩阵R()是上三角矩阵，所以R() 图3 Lyapunov指数谱 R()也是上三角矩阵.因此定义 Fig.3 Lyapunov exponent spectrum (Q()J()(）)g, i>j; 在图3中分别取q=14和9=22得到如图4所示 Kg()= 0, i=j:(12) 的时间序列上的Lyapunov指数曲线. (o()J()0(t)) i<j. 当q=l4时，Lyapunov指数的计算结果如图4(a) 1=4x103 5间 1=-0.06 个 1-0.37 1=-0.11 10 1-6.52 -10 1=-5.32 -15 20 30 40 60 -20 0 20 30 40 50 60 时间s 时间在 图4时间序列上的Lyapunov指数.(a)q=14:(b)q=22 Fig.4 Lyapunov exponent time histories:(a)q=14:(b)q=22工程科学学报，第 38 卷，第 1 期 ( λ ＞ 0) 收缩形成稳定吸引子． 根据 Lyapunov 指数的 性质可知道，只要有一个 Lyapunov 指数为正，就可以 判断该系统处于混沌振荡状态． 系统( 6) 的线性变分方程为 Y · ( t) = J( t) Y( t) ，Y( 0) = I3 ． ( 8) Y∈Ｒ3 × 3，I3是 3 × 3 的单位矩阵，J( t) 是系统的 Jacobi 矩阵 J( t) = f( 1) y1 f( 1) y2 f( 1) y3 f( 2) y1 f( 2) y2 f( 2) y3 f( 3) y1 f( 3) y2 f( 3) y                 3 = 0 1 0 － 1 － κ 1 － β 槡y3 0 － βy1 2 槡y            3  ． ( 9) 由于 Y( 0) 是非奇异的，所以对于 t≥0，Y( t) 也是 非奇异的． Y( t) 可 QＲ 因数分解，记为 Y( t) = Q( t) Ｒ( t) ，上三角矩阵 Ｒ( t) 的对角元素都是正的． 将该 QＲ 分解式带入方程( 8) ，有 Q · ( t) Ｒ( t) + Q( t) Ｒ · ( t) = J( t) Q( t) Ｒ( t) ， Q( 0) Ｒ( 0) = I3 ． ( 10) 将式( 10) 左乘 QT ，右乘 Ｒ － 1，有 QT ( t) Q · ( t) + Ｒ · ( t) Ｒ － 1 ( t) = QT ( t) J( t) Q( t) ， Q( 0) = I3，Ｒ( 0) = I3 ． ( 11) 图 4 时间序列上的 Lyapunov 指数． ( a) q = 14; ( b) q = 22 Fig． 4 Lyapunov exponent time histories: ( a) q = 14; ( b) q = 22 由于 QT ( t) Q( t) = I3，矩阵 K( t) = QT ( t) Q · ( t) 是斜对 角矩阵，并且由于矩阵 Ｒ( t) 是上三角矩阵，所以 Ｒ · ( t) Ｒ － 1 ( t) 也是上三角矩阵． 因此定义 Kij( t) = ( QT ( t) J( t) Q( t) ) ij， i ＞ j; 0， i = j; － ( QT ( t) J( t) Q( t) ) ij． i ＜ { j． ( 12) 从 K( t) = QT ( t) Q · ( t) 中可以得到 Q · ( t) = Q( t) K( t) ，Q( 0) = I3 ． ( 13) 方程( 12) 和方程( 13) 就确定了关于矩阵 Q( t) 的微分 方程． 由于 QT ( t) J( t) Q( t) 是斜对称的，由方程( 11) 得到 Ｒ · '( t) = Ｒ · ii ( t) Ｒii ( t) = ( QT ( t) J( t) Q( t) ) ii，Ｒ' i ( 0) = 0． ( 14) 这里 Ｒ' i ( t) = ln( Ｒii ( t) ) ，Lyapunov 指数的时间演化式 为 λi ( t) = Ｒ' i ( t) /t，所以得到 Lyapunov 指数为 λi = limt→∞ Ｒ' i ( t) t = limt→∞ λi ( t) ，i = 1，2，3． ( 15) 利用 QＲ 分解算法和 MATLAB 软件计算得到如图 3 所示的 Lyapunov 指数谱，λ1 ＞ λ2 ＞ λ3 ． 可以看到， q = 12. 55 时，Lyapunov 指数在形态上发生“突变”，在 章节 2. 1 中分析可知，此时发生擦边分岔． 图 3 Lyapunov 指数谱 Fig． 3 Lyapunov exponent spectrum 在图 3 中分别取 q = 14 和 q = 22 得到如图 4 所示 的时间序列上的 Lyapunov 指数曲线． 当 q = 14 时，Lyapunov 指数的计算结果如图 4( a) · 831 ·