一类基于奇异值分解的Lyapunov指数计算方法

D0I:10.13374/i.issn1001053x.2005.03.029 第27卷第3期 北京科技大学学报 Vol.27 No.3 2005年6月 Journal of University of Science and Technology Beijing Jun.2005 一类基于奇异值分解的yapunov指数计算方法 张晓丹)李志萍》张丽丽) 1)北京科技大学应用科学学院，北京1000832)广东工业大学应用数学学院，广州510090 摘要基于奇异值分解提出了一种有效的计算Lyapunov指数的方法.该方法既适用于离散 系统也适用于连续系统，且避免了频繁的重正交过程.数值计算显示，该方法稳定性好，计算 效率较高，计算结果可靠， 关键词Lyapunov指数：离散系统：连续系统：奇异值分解 分类号0175 对初始值的敏感性是混沌系统的一个最为 X=FX (1) 重要的特征，李雅普诺夫指数就是这种敏感性的 设{X)是该系统的一条轨道，△X为偏离该 量化.1983年，格里波基证明当一系统（非发散） 轨道的一微小量.则△X的演化满足 的李雅普诺夫指数出现正值时，此系统将是混沌 的.此外，判定由混沌系统分解出的子系统间能 Ax-服l小-AX 其中，J为F在X处的Jacobian矩阵，于是，可以得 否实现同步的主要判据之一即子系统的李雅普 出△X4=J△X=JJJ△X.记J=JJ-…J,则 诺夫指数（条件李雅普诺夫指数）皆应为负值，. △X1=△X. 可见计算李雅普诺夫指数谱对于动力系统的混 以下是离散系统Lyapunov指数谱的定义. 沌判定和混沌控制都是极为重要的 定义1对于n维离散动力系统(I),J如上 本文就已知动力学方程的系统，提出了一种 定义.设TX0=1im(W))2,则)是一个正定 有效的基于奇异值分解计算Lyapunov指数的方 矩阵，设≥2…≥>0为T)的n个特征值.则 法，由于该方法是从离散系统出发推得的，故其 系统(l)的第i个Lyapunov指数： 只适用于直接计算离散系统的Lyapunov指数.笔 G,=n2(i=1,2,…,n) (2) 者尝试把连续系统离散化，将其推广到计算连续 1,2连续系统Lyapunov指数的定义 系统的Lyapunov指数.结果表明该方法是完全可 考虑如下n维连续动力系统： 行的.本文详细给出了用定义法计算动力系统 X=F(X) (3) Lyapunov指数的具体实现过程，使Lyapunov指数 其中，X=(x,,x)P∈R,FX)=(f(X),(X,…, 的计算变得简单易行，并举例用这两种方法进行 了数值计算 方约r文=d d 设系统(3)以七为初始条件形成一轨道为 1 yapunov指数的定义 x化，)，若其初始条件的微小变化为(x,0),在1 时刻时轨道的偏差为(x,),则 对于维数大于1的系统，存在着Lyapunov指 Sx(xo,1)=xt,xo+x(xo,0))-x(t,xo). 数的集合，或称为Lyapunov指数谱.每一个Lyap (4) unov指数刻划系统的两条相互靠近的轨道沿某 定义o=lim In,p -=18(,0l 一特定方向按指数相互分离或靠拢的速率， 称o为系统(3)的Lyapunov指数，它表示相邻轨道 l.1离散系统Lyapunov指数的定义 发散（或收敛）率的测度， 考虑如下n维离散动力系统： 下面是连续系统Lyapunov指数谱的定义. 定义2对于n维连续动力系统(3)，在t=0时 收稿日期：2004-05-25楼回日期：200412-27 以x为中心，a(o,0)训为半径作n维球面.由于 基金项目：国家自然科学基金资助项目QNo.70271068) 作者简介：张晓丹(1960一)，女，教授 各方向收缩或扩张程度的不同，随着时间的演

第 2 7 卷 第 3 期 20 5 年 6 月 北 京 科 技 大 学 学 报 JO u rn a l o f U n iv e rs i ty o f s e le n e e a n d 及e h n o l o gy B e ij i n g 、 b L2 7 N o J J U n . 2 0 0 5 一类基于奇异值分解的 Ly aP u n ov 指数计算方法 张 晓丹 ” 李志 萍 ” 张 丽 丽 ” l )北 京科技 大学应 用科 学学 院 , 北京 】0 0 0 83 2) 广 东工 业大 学应用 数学 学院 , 广州 5 10 0 90 摘 要 基 于奇异 值分 解提 出 了一 种有 效的计 算 切叩oun v 指数 的方 法 . 该方 法 既适用 于离散 系统 也适用 于连 续系统 , 且避 免 了频繁 的重正 交过程 . 数值 计算显 示 , 该方法稳 定性好 , 计算 效率 较 高 , 计算 结果 可靠 . 关键 词 yL aP un vo 指 数 ; 离 散系统 : 连续 系统 ; 奇异值 分解 分 类号 0 1 7 5 对 初 始 值 的敏 感 性 是 混沌 系 统 的 一 个最 为 重要 的特 征 , 李雅 普诺 夫指 数就 是这种 敏感 性 的 量化川 . 19 83 年 , 格里波基 证 明当一 系统 ( 非发 散 ) 的李 雅普诺 夫指 数 出现正 值 时 , 此 系统 将是混 沌 的份, . 此 外 , 判定 由混沌 系统 分解 出的子 系统 间能 否 实现 同步 的主 要 判据 之 一 即子 系 统 的 李雅 普 诺 夫指 数 ( 条件 李 雅普诺 夫 指数 ) 皆应 为 负值`圳 . 可 见 计 算 李雅 普 诺 夫指 数 谱对 于 动 力 系统 的混 沌 判 定和 混沌 控 制都 是极 为重 要 的 , 本文 就 已知 动力 学方 程 的系统 , 提 出了一种 有 效 的基 于奇 异值 分 解计 算 yL 叩u n o v 指数 的方 法 . 由于 该方 法 是从 离散 系统 出发推 得 的 , 故其 只 适用 于直接 计 算离 散系 统 的 yL aP uon v 指数 . 笔 者 尝试把 连续 系统 离散 化 , 将其 推广 到计算 连续 系 统 的yL aP u n o v 指数 . 结 果表 明该 方法 是完 全可 行 的 . 本 文 详细 给 出 了用 定 义 法 计算 动 力 系统 yL 即un vo 指数 的具 体 实现 过程 , 使 yL a P un vo 指数 的计 算变 得简 单易行 , 并举例 用这 两种 方法进 行 了数 值 计算 . 戈 辛 , = 月尤 ) ( l ) 设 {Xk } 是该 系统 的一条 轨道 , △戈 为偏 离该 轨 道 的一微 小 量 . 则 △戈 的演化 满足 ~xA 1一 嚼} 、 卜 一 “ 1 yL a P u n o v 指 数 的定义 对 于 维数 大于 1 的系 统 , 存 在着 切 a Pun o v 指 数 的集合 , 或称 为 yL aP un vo 指数 谱 . 每 一个 切叩 - un vo 指数 刻划 系 统 的两 条相 互靠 近 的轨 道沿 某 一 特 定方 向按 指 数相 互分 离 或靠 拢 的速 率 . L l 离散 系统 yL aP u on v 指 数 的定 义 考 虑如 下 n 维 离 散动 力系 统 : 收稿 日期 : 2 0 04 es 习5佗5 修 回 日期 : 2 0 0冬1 2佗 7 基金 项 目 : 国家 自然 科学基 金资助 项 目困.0 7 0 2 71 0 6 8) 作者 简介 : 张 晓丹 ( 1% O一) , 女 , 教授 其 中 , kJ 为F 在龙 处 的 J ac ob i an 矩 阵 . 于 是 , 可 以得 出八阶 毕 : = 工△戈 = 再几 一 , … jJ △X i . 记了” = 去大 一 」 … , 则 八脸 咔 : = 少k) A X i . 以下 是离 散 系统 yL a P u n o v 指数 谱 的定 义 . 定 义 1 对 于 ” 维 离散 动力 系统 ( l) , 产如 上 定义 . 设 爪幻 = l而(幼k’) 丫即2k)K , 则 双幻 是 一个正 定 矩 阵 , 设又 1之儿沙 二 全又户0 为 洲幻 的 n 个特 征值 . 则 系统 ( l) 的第 i 个 yL aP un vo 指数 : 民 = 以 , ( i = l , 2 , … , n ) ( 2) 1.2 连 续 系统 yL a p u n vo 指数 的 定义 考虑 如 下 n 维 连续 动 力系 统 : X = 月幻 ( 3 ) 其 中 , X = (x , , 关 , … , 瓜) T 任 r , 月幻 = 以因 , 关因 , … , ” 了 一 ` 、 , ; , _ d X 几(幻)T , X = 弓于 . J 。 。 二刀 , ` 一 dr ’ 设系 统 (3 ) 以与 为初 始 条 件形 成 一轨 道 为 城t , x0 ) , 若 其初 始条 件 的微 小变 化 为 叙伙。 , 0) , 在 t 时 刻 时轨道 的 偏差 为叙认 。 , )t , 则 叙伙。 , )t = 城t沐`十叙伙` , 0) ) 一城才丙) . 定义 一 吩 场耀黝 称 , 为 系统〔3) 的 yL a PnU vo 指 数 , 发散 ( 或 收敛 ) 率 的测度 . (4 ) 它表 示 相邻 轨道 下 面 是连 续系 统 yL 叩un vo 指数 谱 的定 义 . 定 义 2 对 于 n 维 连续 动力 系统 (3 ) , 在 t = O 时 以 与 为 中心 , }}叙认 。 , 0) }}为半径 作 n 维球 面 . 由于 各方 向收 缩或 扩 张程 度的 不 同 , 随着 时间 的 演 DOI: 10. 13374 /j . issn1001 -053x. 2005. 03. 029

·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 〕酞 全州不 ` , 多 )

Vol.27 No.3 张晓丹等：一类基于奇异值分解的Lyapunov指数计算方法 ·373· 其中，△xr为一微小量.则H在y,处的Jacobian矩阵 2=ind 为+J(x,W)△x,其中，Jx,y)为f在y.处的Jacobian 矩阵.令+x,y)△x=J,则也完全类似于离散动 力系统. =-9+T2n1. 在实际计算中常将d取为1. 3基于定义2的算法实现（方法二） 4 数值计算 具体算法为：将(，0川取为d(d为常数)， 以为球心，欧几里德范数为d的正交矢量集{e, (I)Henon映射. x1=1-1.4x+y r,,e,}为初始球.由非线性微分方程(3)分别计 y:=a3 算出点xo,+e,+e,…,+e,经过时间T演化的轨 其Lyapunov指数的理论值为0.418，-1.622， 迹，设其终了点分别为，oi,…,oa.令"=xo- 用“方法一”计算其Lyapunov指数，结果为 o,8=2一x,…,8=m一x0,则得到新的矢量 0.4162,-1.6201,表明该方法是完全可行的. 集{8r",8x,…,8c}. 下面分别用“方法一”、“方法二”和QR分解 由于各矢量在演化过程中都向最大Lyapun- 法计算Lorenz系统、Chen系统和Rossler系统的 ov指数方向靠拢，所以需要通过Schmidt正交化 全部Lyapunov指数.当然，用“方法一”和QR分 不断对新得到矢量集进行置换.即 解法计算连续系统的Lyapunov指数需先将其离 =ax, 散化. 的 (2)Lorenz系统. =6"-(6,, [i=a(y-x) 位=cx-y-xz 4 =xy-bz 当=16，b=4,c=45.92时，该系统呈现混沌性 =8-(6x","a"-(a",》"-- 态.表1是用“方法一”、“方法二”和QR分解法计 (",", 算该系统Lyapunov指数的结果比较. 4。网 (3)Chen系统. [i=a(y-x) 接着以m为新球心，范数为d的正交矢量集 =(c-a)x-xz+cy {du",d,,du"}为新球继续上述演化.设演化 2=-bz+xy 到第N步时得到矢量集{”，”，，m,并且N足 当a=35,b=3,c=28时，该系统呈现混沌性态. 够大，则可以得到连续系统(3)Lyapunov指数的近 表2是用“方法一”、“方法二”和Q分解法计算 似计算公式： 该系统Lyapunov指数的结果比较. 表1计算Lorenz系统yapunov指数的结果比较 Table 1 Comparison of Lyapunov exponents in the Lorenz system 理论值 方法一 QR分解法 方法二 4阶R-K离散 欧拉离散 4阶R-K离散 欧拉离散 1.5376.-0.0074, 1.5124,0.0216 1.4903,0.0325. 1.5135,-0.0437, 1.4950,-0.1115 1.50,0,-22.5 -22.5301 -22.5944 -22.5236 -22.4698 -22.4439 误差% 0.22 0.43 0.18 0.24 0.55 表2计算Chen系统yapnnov指数的结果比较 Table 2 Comparison of Lyapunov exponents in the Chen system 方法一 QR分解法 理论值 方法二 4阶R-K离散 欧拉离散 4阶R-K离散 欧拉离散 2.0602,-0.0129, 1.8078,0.0075, 1.9398.-0.1449, 2.0579.-0.0159 1.7448.0.0381. 2.0274,0,-12.0135 -12.0731 -11.6565 -11.7947 -12.0420 -11.6391 误差% 0.57 3.44 2.27 0.37 3.86

认〕1 . 2 7 N o . 3 张 晓丹等 : 一类 基于 奇异值 分解 的 Ly aP un ov 指 数计 算方法 . 3 73 . 其 中 , xA 为一 微 小量 . 则 H 在y `处 的 J ac ob i an 矩 阵 为介浏禹小)公 , 其 中 , 洲元小 )为 f 在 y ` 处 的 J ac o ih an 矩 阵 . 令 介j 吮 亩, yt) 酝二 试 , 则 也完全 类似 于 离 散动 力系 统 . 3 基 于 定义 2 的算法 实 现 (方 法 二 ) 具体 算 法 为 : 将 }}叙编 , 0)l }取 为 d d( 为 常数 ) , 以x0 为球 心 , 欧几 里德 范数为 d 的正 交矢 量集 { el, 几 , … , en } 为初 始球 . 由非 线 性微 分方 程 (3) 分别 计 算 出点与 , x0 十 el , x0 十氏 , … , 与+en 经过 时间 厂演化 的轨 迹 , 设其 终 了点分 别 为 xo , x0 : , … , x0 。 . 令酬 , , 二 x0 1一 xo , 碳 l) 二 x0 2 一 xo , … , 酬 〕 =x0n 一 x0 , 则得 到 新 的矢 量 集 {酬 ” , 酬 1) , … , 聪 , ) } . 由于各 矢量 在 演化 过程 中都 向最 大 yL aP un - vo 指 数 方 向靠拢 , 所 以需 要通 过 S e知rn i dt 正 交 化 不 断对 新得 到 矢量 集进 行 置换 . 即 城 , , = 叙} , , , 、 , 一半 + 命誉vln4I 毗 、 。 一粤 十命誊 , n ,,诸) !, · 在 实 际计 算 中常 将 d取 为 4 数值计 算 ( l ) He n o n 映射 . 1一 1一 1一 1 . 叫+y, 口 X3 , 试 , , “ , ’ - 下(不 ’ 谓, = 酬 ,一 (议 1) , u } , , ) u } , , , ;vl , u , = 不可 , 其 yL ap un vo 指 数 的理 论值 为 .0 4 18 , 一 1 . 6 2 , 用 “ 方 法 一 ” 计 算 其 yL ap un vo 指 数 , 结 果 为 .0 4 1 62 , 一 1 . 62 0 1 , 表 明该方 法是 完 全可 行 的 . 下面 分别用 “ 方 法 一 ” 、 “ 方 法 二 ” 和 QR 分解 法「幻计 算 L o er nz 系 统 、 C h en 系 统和 R 6 s ler 系统 的 全 部 yL ap un vo 指 数 . 当然 , 用 “ 方 法一 ” 和 QR 分 解 法计 算 连续 系 统 的 yL aP un vo 指 数 需先将 其 离 散 化 . (2 ) L or enz 系统 . }名 } 几万Z = a 伽一 x ) = 〔沈 一〕 夕一戈 2 = 划一 bz 谓 ,= 酬 ,一 (衅 ’ ,斌 , , (以i ) , u 华 , > u 华 : , >衅,一 (麟 l ) ,对〕 )砖,一 一 _ “ , 刃 一 刀可下 . 当二 16 , b科 , 。 = 4 5 . 9 2 时 , 该 系统 呈现 混沌 性 态 . 表 l 是 用 “ 方 法一 ” 、 “ 方法 二 ” 和 QR 分解 法计 算 该 系统 yL ap un vo 指 数 的结 果 比 较 . ( 3 ) C h e n 系 统 . 接 着 以xo 为 新 球 心 , 范 数 为 d 的 正 交 矢 量 集 { d衅 ,d 谓 ) , … , d川, }为新 球 继续 上述 演 化 . 设演 化 到第 N步 时 得到 矢量 集 {砂 , 谬 , 一 , 刃 } , 并且 N 足 够大 , 则 可 以得 到连 续 系统 ( 3) yL aP ~ v 指 数 的近 似 计 算 公式 : }舍 } t省Z = a 伽一 x ) = ( c 一 a )x 一脂+ cy = 一 西之十郑夕 当a 二3 5 , b二 3 , 份2 8 时 , 该 系统 呈现 混 沌性 态 . 表 2 是用 “ 方法 一 ” 、 “ 方 法二 ” 和 QR 分 解法 计算 该 系统 yL ap un vo 指 数 的结 果 比 较 , 表 1 计 算 L O代nz 系 统 yL aP un vo 指 数 的结果 比 较 aT bl e 1 C o m P a ir , o n o f yL a P u n vo e x op n e n st i . 比e L o er n z s ys t e m 理 论值 方 法一 QR 分解法 4 阶 R es K 离散 欧拉离 散 方法二 1 . 50 , 0 , 一 2 2 . 5 1 . 5 3 7 6 , 一 0 . 0 0 7 4 , 一 2 2 . 53 0 1 1 5 12 4 , 0 . 0 2 1 6 , 一 2 2 . 59 4 4 1 . 49 0 3 , 0 . 0 3 2 5 , 一 2 2 5 2 3 6 4 阶 卜K 离散 1 . 5 13 5 , 一 0 . 0 4 3 7 , 欧拉 离散 一 22 4 6 9 8 1 . 4 9 5 0 , 一 0 , 1 11 5 一 2 2 4 4 3 9 误差 /% 表 2 计算 C b e n 系 统 yL a P u o o v 指 数的 结果 比较 aT b l e 2 C o m P a d s o . o f yL a P u n o v e x P o n e n tS i n ht e C b e n s ys t e m 理 论值 方法 一 方 法二 Q R 分解 法 4 阶 R es K 离 散 欧 拉离散 2 乃2 7 4 , 0 , 一 1 2 0 13 5 2 . 0 60 2 , 一 0 0 12 9 , 一 12 刃7 3 1 1 . 80 7 8 , 0 0 0 7 5 , 一 11 . 6 5 6 5 1 . 9 3 9 8 , 一 0 . 14 4 9 , 一 1 1 . 79 4 7 4 阶 R , K 离 散 2 . 0 5 7 9 , 一 0刃 15 9 欧拉离 散 一 12 . 04 2 0 1 74 4 8 , 0 . 03 8 1 , 一 1 1 63 9 1 误 差 肠

374· 北京科技大学学报 2005年第3期 (4)Rossler系统， 当a-0.15,b-020,c-10.0时，该系统呈现混沌 x=-y-2 性态，表3是用上面三种方法计算该系统Lyapun- y=xtay ov指数的结果比较， =b+xz-cZ 表3三种方法计算Rossler系统Lyopunov指数的结果比较 Table 3 Comparison of Lyapunov exponents in the Rossler system calculated by the above three methods 方法一 QR分解法 理论值 方法二 4阶R-K离散 欧拉离散 4阶R-K离散 欧拉离散 0.0903,0.0383 0.0937.0.0408 0.0865,0.0221 0.0861,0.0423 0.0895.0.0448 0.09.0,-9.79 -9.7503 -9.8057 -9.7372 -9.7501 -9.8055 误差% 0.56 0.45 0.59 0.60 0.48 5 结论 Phys Rey Lett,1991,A44:2374 [5)盛昭翰，马军海.非线性动力系统分析引论。北京：科学 基于奇异值分解的计算Lyapunov指数的方 出版社，2001 [6)王东生，曹磊。混沌分形及其应用，合肥：中国科学技术 法计算连续系统在多数情况下略优于QR分解 大学出版社，1995 法：与定义法相比较，避免了频繁的重正交过程， [刀刘秉正.非线性动力学与混沌基础.长春：东北师范大学 且计算结果靠近于理论值（不同的时间步长会使 出版社，1995 计算结果有一定偏差).在求离散系统yapunov [8]Hubertus F.von Bremen,et al.An efficient QR based method for the computation of Lyapunov exponents.Physica,1997,D101:1 指数时，步骤简单且计算精度较高.总之，该方法 [9]Zhang X D,Min L Q.Theory for constructing generalized syn- 不失为一种有效的计算动力系统Lyapunov指数 chronization and application.J Univ Sci Technol Beijing,2000, 的方法. 7(3):235 [10]Zhang X D,Min L Q.New theorems on generalized synchron- 参考文献 ization of chaos and applications.Adv Syst Sei Appl,2000, [1]Heinx O P.Chaos and Fractals.Newyork:Springer Verlag,1992 (2:37 [2]John A.An Exploration of Chaos.North-Holland:Amsterdam, (11]Zhang X D,Zhang LL,Min L Q.Construction of generalized 1994 synchronization for a kind of array differential equations and ap- [3]Pecore L M,Carroll T L.Synchronization in chaotic system plications.Chin Phys Lett,2003,20(12):2114 Phys Rev Lett,1990,64:821 12]张晓丹，张丽丽，闵乐泉.阵列动力系统广义同步的若干 [4]Pecore L M,Carroll T L.Driving systems with chaotic signals. 新理论及计算机模拟.北京科技大学学报，2004,26(2)：211 A method based on singular value decomposition for computation of Lyapunov ex- ponent ZHANG Xiaodan,LI Zhiping",ZHANG Lili 1)Applied Science School.University of Science and Technology Beijing.Beijing 0008.China 2)School of Applied Mathematics,Guangdong University of Technology,Guangzhou 510090,China ABSTRACT An efficient method based on Singular Value Decomposition(SVD)for the computation of Lyapunov exponent is presented.This method is applicable not only to discrete systems but also to continuous systems,and it does not require frequent reorthogonalization.Numerical calculations show that the method is stable and has a hig- her computing efficiency,and the results are credible. KEY WORDS Lyapunov exponent,discrete system,continuous systems,singular value decomposition

3 7 4 . ( 4 ) 助 s s l e r 系 统 . )t Z 二 一夕一z 二 x + ay 二 b 斗xz 一 cz 北 京 科 技 大 学 学 报 2 0 5 年 第 3 期 当a 二0 . 1 5 , b = .0 2 0 , 二 1.0 0 时 , 该 系统呈 现 混沌 性 态 . 表 3 是 用上 面三 种方 法计 算 该系统 yL ap un - vo 指 数 的结 果 比 较 . 表 3 三 种方 法计算 R 6s l盯 系 统 yL 叩un vo 指数 的结果 比较 aT bj e 3 C o m P a ir s o n o f yL a p u n vo e x P o . 叨 t s i n 比e 助 s s l e r , y s t e m e a k u l a tde b y t h e a b vo e t h比 e m et 七o d s 方法 一 理 论值 4 阶 体K 离 散 0乃9 0 3 , 0 . 0 3 8 3 欧拉 离散 方法二 Q R 分解 法 0 . 0 9 . 0一 9 , 7 9 0 . 0 9 3 7 , 0 乃40 8 0 乃86 5 , 0刀2 2 1 4 阶 几K 离散 0刀8 6 1 , 0 0 4 2 3 欧拉离 散 0 刀8 9 5 , 0 , 0 4 4 8 一 9 . 7 50 3 一 9 名0 5 7 一 9 夕3 7 2 一 9 . 7 50 1 一 9 名0 5 5 误 差 /% 0万6 0 4 5 0 乃9 0 . 60 0 .4 8 5 结 论 基 于奇 异值 分解 的计 算 yL aP un ov 指 数 的方 法计 算 连续 系统 在 多数 情况 下 略优 于 QR 分解 法 ; 与定义法 相 比 较 , 避 免 了频繁 的重正 交过程 , 且 计 算结 果靠 近于 理论值 ( 不 同的 时间步 长会使 计算 结果 有一 定偏 差 ) . 在 求离 散系 统 yL aP un vo 指数 时 , 步骤 简单 且计算 精度较 高 . 总之 , 该方 法 不 失为一 种有 效 的计算 动 力系 统 切aP un vo 指 数 的方 法 . 参 考 文 献 [ 11 H e i n x 0 P C h ao s an d F r ac at l s . N e wy o kr : S P inr g e r Ve r lag , 19 9 2 [ 2」J o hn A . nA EXP l o r at i o n o f C h ao s . N o hrt 一 H o llan d : A m st e dr am , 19 9 4 [ 3 ] P e e o r e L M , C ar o ll T L . S y n e h r o n i atZ i o n i n e h ao t i e sy s t e m . P h y s R vc L e t , 19 9 0 , 64 : 8 2 1 1 4 」 P e e ore L M , C a了r o ll T L . D r i v ing sy ste m s w iht e h aot i e s ig n a l s - P b y s eR v L e t , 19 9 1 , A 44 : 2 3 7 4 【5] 盛昭 翰 , 马军海 . 非线性 动力系 统分析 引论 . 北京 : 科学 出版社 , 2 00 1 6[ ] 王 东生 , 曹 磊 . 混 沌分 形及其应 用 . 合 肥 : 中国科 学技术 大学 出版社 , 19 9 5 7[] 刘秉正 . 非线性 动力学 与混沌 基础 . 长春 : 东北师 范大学 出版社 , 19 9 5 [ 8】H u b e rtU s F vo n B re m e n , et al . nA e if e i e in Q R bas e d m het o d ofr hte e o m P u t a t ion o f yL 叩un o v e XP o n e nts . P h y s i e a , 19 9 7 , D 10 1 : l [9 1 Z h an g X D , M in L Q . hT e o ry ofr e o ns trU c tin g g e n e ar lize d s yn - e h r o n i atZ i o n an d aP P li e at i o n . J U n 扮 S e i eT c h . o l B e ij i . g , 2 0 0 , 7 ( 3) : 2 3 5 [ I 0 1 Z h an g X D , Mi n L Q . N e w ht e o re m s o n ge n era l i z e d s yn c h r o n - i atZ i o n o f e h a o s an d ap Pli c a ti on s . A d y S ys t S e i A p p l , 2 0 0 , ( 2) : 3 7 f l l l Z h an g X D , Z h an g L L , M i n L Q . C o n s如 e ti佣 o f g e n aer lize d s y n e h r o n i azt i o n ofr a k in d o f 日r a y d i fe 溯 t ial qe u a tion s an d 即 - P li e at i o n s . C h i n 外羚 L e t , 2 0 0 3 , Z D` 12 ) : 2 114 〔12 ] 张晓丹 , 张 丽丽 , 阂 乐 泉 . 阵列动力 系统广 义同步 的若干 新理论及计算 机模拟 . 北京科技 大学学报 , 2 0 04 , 26 ( 2) : 21 1 A m e t h o d b a s e d o n s i n g u lar v a l u e d e e o m P o s i t i o n of r e o m P u t a ti o n o f yL aP un o v e x - P o fl e l t 乙从咬刃 G 方飞a o da 月 l ’ , lL hZ iP 矛双9 1’ , 刀丈扭刃G L ilz , 1) A p li e d S e i e n e e S e h o o l , U n i v ers ity o f s e i e n e e an d eT e hn o l o 留 B e ij i n g , B e ij ing l 0 0 0 83 , C h in a 2 ) S e h o l o f A P Pli e d M hat em at i e s , G u an g d o gn U n i v e rs ity o f 介 e hn o l o gy, uG an ghz ou 5 I0 0 9 0 , C h ina A B S T R A C T nA e if e i e in m e ht o d b as e d o n S i n g u lar Va l u e D e e o m P o s it i o n (S V D ) fo r ht e e o m P ut at i o n o f yL aP uno v e xP o n e nt 1 5 rP e s e nt e d , hT i s m e ht o d 1 5 a P P li e ab l e n o t o n ly t o d i s er e te sy s t em s b ut al s o to c o nt l n u o u s s y st e m s , an d i t d o e s n o t er qiu r e fr e q u e in r e o rt h og o n a li atZ i o n . N um e ir e al e a l cul at ion s s h o w ht at ht e m e ht o d 1 5 s tab l e a n d h as a hi g - he r e o m P u tign e if e i e n e y, an d th e er s u lt s aer e er d i b l e . K E Y W O R D S yL aP un o v e x P o n e n t , d i s e er te s y s et m , e on t in uo us s y s et m s , s i n g u lar v a l u e d e e o m P o s i t i o n