·372· 北京科技大学学报 2005年第3期 化，在1时刻此球面将变形为n维椭球面.设此椭 1,2,",n.于是系统(1)的Lyapunov指数为o= 球面第i个坐标轴方向的半轴长为8x,,训，则 imnn=lims亿+ina,l,2n.当k较大 系统(3)的第i个李雅普诺夫指数a,为 n1nl8x(o,t训 时，a≈s空na,=1,2n.具体算法为： a,=lim nx( (i=1,2,…,n (5) 。=5S=lnm 2定义1的算法实现（方法一） b=0,e=(1,1,…,1)r For i=1 to k 定理1设J=JJ,…J,若J“的奇异值为 B,=J.V-IS- i=1,2,…,n),则系统(1)的Lyapunov指数为o,= [V S,U.]=svd(B.), mnn=1,2,. a,=max(diag(S)), 证明设的奇异值分解式为=VS:U, 5=d 其中，V,U为正交矩阵，S=diag,2,,n,则 b=b+log(a.) ()=(U)SU.令T:=(JT)2,则T= End(For循环结束) (U)SU.即与S相似，有相同的特征值.而 LCEvector=[log(diag(S))+bel/k S=diag(i,,,).故T的特征值为(i=1,2,", 对于连续动力系统，若用此方法求其Lyapun- n),从而T的特征值为n(i=1,2,…,n). OV指数，可以先进行离散化，将其转化为离散动 由定义1，T)=limT,T八)的特征值为λ，= 力系统.当然，将连续系统离散化的方法很多.这 lim"(=l1,2,,).从而系统(I)的yapunov指数 里，给出两种离散化的方法：四阶Runge-Kuta法 和欧拉法. a=M,=lmn.(i=l,2,,m小.基于定理1对 设有一向量微分方程为： 进行奇异值分解，令V,=S=l,I为n阶单位矩阵). y=f(x,y) svd(J)=svd[J.J-…J]=svd[JJ-…JJVS)]= 其中，x∈R,y=(W",y,…,y,fxy=[f"(x,y以， svdJ-…J(JVS)]U]=… svd[JJ-…(J,S-j[Un-…U]=… x以，xjj= dx VS[U…UVS,U. (a4阶R-K方法离散.采用四阶R-K方法可 其中，V,U,U均为正交阵，S为对角矩阵，VSU 将其离散为如下形式： 为J'S-:的奇异值分解，的奇异值即为S的 y-y.+G(K+2K.+2K+K)). 主对角线元素. 其中，K=x小，K=x+资+hK. 由于在上述分解过程中出现J与S(=1,2,, K-f(x.+2.y.+hK).K.=/(x.th.y.+hk). k-1)的累次相乘，而且实际计算中k又相当大， 可能导致计算滋出，故分解过程中采用如下修正 令=Jx小，h=Jx+空，以+hK 方法：对J"-S-i=1,2,,k-1,S=D进行奇异值 万=x+空y+山=九hK. 分解，JYS=VS,U,然后将S,中绝对值最大的 其中，J为f对y的Jacobian矩阵，G在y,处的Jaco 元素提出（实际上，S,(对角矩阵）中主对角线上的 bian矩阵为： 元素均为正值，故只要将其最大元素提出即可)， 1+gJ+22J+242hW+2Wj+ 从而使得提出后新得到的矩阵5，中的每一元素 都不大于1，即Jy-S.-=a,S,U.(其中a,为S的最 J+hJ[+h+2J0]》. 大元素).则此时的分解过程为： 其中，I为单位阵，计算时只要令 svd)=svd[J.J-…(JVSa)】]= 1+gU+2J0+2MH2J+2W42W+ a]svdJJJ(VS)U= J.(I+hJ[I+]hJ:(I+]hJ)])}=J.. {aa]svd[J:J.-…(JV52)][UU1]=… 则完全类似于离散动力系统， [aa…ak-lVSU…U](Vo=S=0. (b)欧拉法离散.采用欧拉法可将其离散为如 其中，a为S,0=1,2,,k-1)的最大元素.由上面 下形式： 式子可得到J的奇异值为n=aaa-S(i,),户 -y+fx,y)△H(xy,）.北 京 科 技 大 学 学 报 2 0 0 5 年 第3 期 化 , 在 t 时刻此 球面将 变形 为 n 维椭球 面 . 设此椭 球面 第 i 个坐 标轴 方 向的半轴 长为 }{阮俩 , )ltI , 则 系统 ( 3) 的 第 i 个 李雅 普 诺夫 指数 ia 为 2 , … , n . 于 是 系 统 ( l) 的 yL a P u n vo 指 数 为厂 所、 卿和*nS i(, i) + 艺 I n马 ) , =1 1 , 2 , … , n . 当k较大 氏 一 ! : 令 ` n }黯黝 ( ,一 ` , 2 , 4 ~ , · , ( 5) 时 , 。 二 冬l( skn i(, i) + 窗nI 。 ) , 、 1 , 2 , … , , . 具体 算 法为 : 凡 产 1 2 定 义 1 的算法 实 现 ( 方法 一 ) 定 理 1 设尹 , = 大大 一 1 … 试 , 若 j ” 的奇 异值 为 袱i=l , 2 , … , n) , 则 系统 ( 1) 的 yL aP un vo 指 数为氏 = 卿士 , n。 ( `一 1 , 2 , … , 。 ) . 证 明 设尹 , 的奇 异 值 分解 式 为尹 , “ 从及 队 , 其 中 , K , 队 为 正 交 矩 阵 , S * = di眼切 , , 叮2 , … , 峥 。 ) , 则 (j ` , )丫 ` , = (从 )议认 . 令 爪= (夕 ` , )丫 ` , ) ” `, ` , , 则 暇 k = ( 队)叹 队 . 即 双 `与瑟相似 , 有相 同 的特 征 值 . 而 瑟 二id ag (呱书 , … , 祠 . 故 鲜 ` 的特 征值 为树( i=l , 2 , … , n ) , 从而 爪 的特 征 值 为叮刀 ` ( i = l , 2 , … , n ) . 由 定 义 1 , 八幻 = h m 工 , 联幻 的特 征 值 为又二 li m 叮) z` ( i = l , 2 , … , n ) . 从 而 系 统 ( l ) 的 琢即un o v 指 数 走~ 沉 。 一 1川 一 卿十 in 。 ( ,一 1 , 2 , … , n) . 基 于 定理 l 对 J k) 进 行奇异 值分 解 , 令 从 = 国斌 (I 为 n 阶 单位 矩 阵) . vs d( J 与= vs d 〔大大 一 ! … J :」= vs d 〔人kJ 一 : … 远(禹叭凡) ] 二 vs d 〔入大 一 、 …禹(zJ 卜况)〕〔鱿卜 … vs d[ 人人 一 ! … (试厂 一 1况 一 , )〕〔鱿 一 1… 鱿」= … 卜议乳「以… 鱿〕全碑么从 kU . 其 中 , K , 鱿 , 队 均 为 正 交阵 , 戈为对角 矩 阵 J 飞又认 为丈从 一 1瑟 一 : 的奇异 值分 解 . kJ ,的奇异 值 即 为二的 主 对 角线 元素 . 由于 在上 述 分解 过程 中 出现冻 ,与况(=1 1 , 2 , … , k 一 1) 的累次相 乘 , 而且 实 际计算 中 k 又相 当大 , 可 能导致计 算溢 出 , 故分解 过程 中采用 如下修 正 方 法 : 对试阵 1民 一 1 i( = 1 , 2 , … , k 一 1 , s0 = 刀进行 奇异 值 分解 , 试E 一况 一 1二 卜万 C , 然 后将忍中绝 对值 最大 的 元 素提 出( 实际 上 , S ( 对角 矩阵 ) 中主 对 角线 上 的 元 素均 为正 值 , 故只 要 将其 最大 元素 提 出即 可 ) , 从 而 使得 提 出后 新 得 到 的矩 阵民中 的每 一 元 素 都 不大 于 1 , 即试卜 一 l民 一 , = a J 侣 万双 (其 中 a 为况的最 大元 素 ) . 则 此 时的分 解 过程 为 : vs d (J 勺= vs d〔人大 一 、 … (lJ 1 2 吉 ) )卜 「a l 〕 s v d! kJ kj , … 丈(丈叱5 1 )〕「U ! ] = 〔 a l a 2」s v d〔人丈 一 1 … (工矶及) ][乱 ` 1卜… at[ 处 … 氏 一 , I ,修 走 〔以 … u l ] (v0 一 s0 胡 . 其 中 , 马 为况心二 1 , 2 , … , k一 l) 的 最大 元素 . 由上 面 式 子 可 得 到尹 , 的 奇 异 值 为粉于al a Z… 久 一 1凡i(, i) , 卜 践二 s0 城 . x 。 b = 0 , e = ( l , l , … , l ) T F or =1 1 ot k B , 二 试E _ .况 一 : , 【E , 民 , 以] = s v d (B , ) , a , = m a x (d iag (民)) , 瓦一 气“ - b = b + 1 0 9 ( a , ) En d ( F o r 循 环 结束 ) L C E v e e t o 二〔1 0 9 ( di a g (从)) + b e 」k/ 对 于连续 动 力系统 , 若用 此方法 求其 yL aP un - o v 指数 , 可 以先进 行 离散 化 , 将其 转化 为 离散 动 力系统 . 当然 , 将 连续 系统离 散化 的方法 很多 . 这 里 , 给 出两 种离 散化 的方 法 : 四阶 R u n g e . K u t ta 法 和 欧拉法 . 设 有一 向量 微分 方程 为 : y 气八才 ,夕) 其 中 , x E R , 少= 甲 , , , 少 , , , … ,夕 ” , ) T , . 八无少) = 了 ” (x , 少) , 厂 ” (x ` , , 一尸x( , ,〕 T , 夕 一 会 . (a) 4 阶 R , K 方法 离散 . 采 用 四阶 Res K 方法 可 将其 离散 为如 下 形式 : , 附 l 一 , 。 如 +二 + 2二二)。 、 , 。 ) . 其 中 , ` =fx(n ,二 , , 二孤 , 谙 , 苏十合斌 ) , 二 =f(x 岭 , yn小z)K , 二 一 f(xn +h , , 汁 hsK ) . 令` 一 x(J ” 二 , , ` 一 (xnJ 垮 , , 。 号斌 . ) , 姚 一 (xJ n 垮 , 苏号hKz ) , ; 一 x("J + h, 汁 h,K ) . 其 中 , J 为f 对 y 的 Jac ob ian 矩 阵 , ` 在苏处 的 J ac o - ib an 矩 阵 为 : 暗 {.l + 2丈仍合M ) + 二。叶 。*叶、 ): + ; +{I hjs t叶hzJ 叶hlJ )}] . 其 中 , I 为 单位 阵 . 计算 时只 要 令 蜡+{J, 2zJ (,4 ,hJ +) s2J 〔叶hzJ叶hlJ )+l 人+{, haJ 卜告从叶h,J )〕}卜 ; , 则完 全类 似于 离散 动 力系统 . (b ) 欧拉法 离散 . 采用 欧拉法 可将其 离散 为如 下 形 式 : y o l = y次f’( x i, y 〕酞 全州不 ` , 多 )