matlab求最大李雅普诺夫(Lyapunov)指数程序
-Lyapunov指数的计算方法
【总结】Lyapunov指数的计算方法非线性理论近期为了把计算LE的一些问题弄清楚，看了有7～9本书！下面以吕金虎《混沌时间序列分析及其应用》、马军海《复杂非线性系统的重构技术》为主线，把目前已有的LE计算方法做一个汇总！1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。
关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍内容。
（1）定义法定义法求解Lyapunov指数.JPG关于定义法求解的程序，和matlab板块的“连续系统LE求解程序”差不多。
以Rossler系统为例Rossler系统微分方程定义程序function dX = Rossler_ly(t,X)% Rossler吸引子，用来计算Lyapunov指数% a=0.15,b=0.20,c=10.0% dx/dt = -y-z,% dy/dt = x+ay,% dz/dt = b+z(x-c),a = 0.15;b = 0.20;c = 10.0;x=X(1); y=X(2); z=X(3);% Y的三个列向量为相互正交的单位向量Y = [X(4), X(7), X(10);X(5), X(8), X(11);X(6), X(9), X(12)];% 输出向量的初始化，必不可少dX = zeros(12,1);% Rossler吸引子dX(1) = -y-z;dX(2) = x+a*y;dX(3) = b+z*(x-c);% Rossler吸引子的Jacobi矩阵Jaco = [0 -1 -1;1 a 0;z 0 x-c];dX(4:12) = Jaco*Y;求解LE代码：% 计算Rossler吸引子的Lyapunov指数clear;yinit = [1,1,1];orthmatrix = [1 0 0;0 1 0;0 0 1];a = 0.15;b = 0.20;c = 10.0;y = zeros(12,1);% 初始化输入y(1:3) = yinit;y(4:12) = orthmatrix;tstart = 0; % 时间初始值tstep = 1e-3; % 时间步长wholetimes = 1e5; % 总的循环次数steps = 10; % 每次演化的步数iteratetimes = wholetimes/steps; % 演化的次数mod = zeros(3,1);lp = zeros(3,1);% 初始化三个Lyapunov指数Lyapunov1 = zeros(iteratetimes,1); Lyapunov2 = zeros(iteratetimes,1); Lyapunov3 = zeros(iteratetimes,1);for i=1:iteratetimestspan = tstart:tstep:(tstart + tstep*steps); [T,Y] = ode45('Rossler_ly', tspan, y);% 取积分得到的最后一个时刻的值y = Y(size(Y,1),:);% 重新定义起始时刻tstart = tstart + tstep*steps;y0 = [y(4) y(7) y(10);y(5) y(8) y(11);y(6) y(9) y(12)];%正交化y0 = ThreeGS(y0);% 取三个向量的模mod(1) = sqrt(y0(:,1)'*y0(:,1));mod(2) = sqrt(y0(:,2)'*y0(:,2));mod(3) = sqrt(y0(:,3)'*y0(:,3));y0(:,1) = y0(:,1)/mod(1);y0(:,2) = y0(:,2)/mod(2);y0(:,3) = y0(:,3)/mod(3);lp = lp+log(abs(mod));%三个Lyapunov指数Lyapunov1(i) = lp(1)/(tstart);Lyapunov2(i) = lp(2)/(tstart);Lyapunov3(i) = lp(3)/(tstart);y(4:12) = y0';end% 作Lyapunov指数谱图i = 1:iteratetimes;plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)程序中用到的ThreeGS程序如下：%G-S正交化function A = ThreeGS(V) % V 为3*3向量v1 = V(:,1);v2 = V(:,2);v3 = V(:,3);a1 = zeros(3,1);a2 = zeros(3,1);a3 = zeros(3,1);a1 = v1;a2 = v2-((a1'*v2)/(a1'*a1))*a1;a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;A = [a1,a2,a3];计算得到的Rossler系统的LE为———— 0.063231 0.092635 -9.8924Wolf文章中计算得到的Rossler系统的LE为————0.09 0 -9.77需要注意的是——定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象。
-Lyapunov指数的计算方法
【总结】Lyapunov指数的计算方法非线性理论近期为了把计算LE的一些问题弄清楚，看了有7～9本书！下面以吕金虎《混沌时间序列分析及其应用》、马军海《复杂非线性系统的重构技术》为主线，把目前已有的LE计算方法做一个汇总！1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。
关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍内容。
（1）定义法定义法求解Lyapunov指数.JPG关于定义法求解的程序，和matlab板块的“连续系统LE求解程序”差不多。
以Rossler系统为例Rossler系统微分方程定义程序function dX = Rossler_ly(t,X)%Rossler吸引子，用来计算Lyapunov指数%a=0.15,b=0.20,c=10.0%dx/dt = -y-z,%dy/dt = x+ay,%dz/dt = b+z(x-c),a = 0.15;b = 0.20;c = 10.0;x=X(1); y=X(2); z=X(3);% Y的三个列向量为相互正交的单位向量Y = [X(4), X(7), X(10);X(5), X(8), X(11);X(6), X(9), X(12)];% 输出向量的初始化，必不可少dX = zeros(12,1);% Rossler吸引子dX(1) = -y-z;dX(2) = x+a*y;dX(3) = b+z*(x-c);% Rossler吸引子的Jacobi矩阵Jaco = [0 -1 -1;1 a 0;z 0x-c];dX(4:12) = Jaco*Y;求解LE代码：% 计算Rossler吸引子的Lyapunov指数clear;yinit = [1,1,1];orthmatrix = [1 0 0;0 1 0;0 0 1];a = 0.15;b = 0.20;c = 10.0;y = zeros(12,1);% 初始化输入y(1:3) = yinit;y(4:12) = orthmatrix;tstart = 0; % 时间初始值tstep = 1e-3; % 时间步长wholetimes = 1e5; % 总的循环次数steps = 10; % 每次演化的步数iteratetimes = wholetimes/steps; % 演化的次数mod = zeros(3,1);lp = zeros(3,1);% 初始化三个Lyapunov指数Lyapunov1 = zeros(iteratetimes,1);Lyapunov2 = zeros(iteratetimes,1);Lyapunov3 = zeros(iteratetimes,1);for i=1:iteratetimestspan = tstart:tstep:(tstart + tstep*steps); [T,Y] = ode45('Rossler_ly', tspan, y);% 取积分得到的最后一个时刻的值y = Y(size(Y,1),:);% 重新定义起始时刻tstart = tstart + tstep*steps;y0 = [y(4) y(7) y(10);y(5) y(8) y(11);y(6) y(9) y(12)];%正交化y0 = ThreeGS(y0);% 取三个向量的模mod(1) = sqrt(y0(:,1)'*y0(:,1));mod(2) = sqrt(y0(:,2)'*y0(:,2));mod(3) = sqrt(y0(:,3)'*y0(:,3));y0(:,1) = y0(:,1)/mod(1);y0(:,2) = y0(:,2)/mod(2);y0(:,3) = y0(:,3)/mod(3);lp = lp+log(abs(mod));%三个Lyapunov指数Lyapunov1(i) = lp(1)/(tstart);Lyapunov2(i) = lp(2)/(tstart);Lyapunov3(i) = lp(3)/(tstart);y(4:12) = y0';end% 作Lyapunov指数谱图i = 1:iteratetimes;plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)程序中用到的ThreeGS程序如下：%G-S正交化function A = ThreeGS(V)% V 为3*3向量v1 = V(:,1);v2 = V(:,2);v3 = V(:,3);a1 = zeros(3,1);a2 = zeros(3,1);a3 = zeros(3,1);a1 = v1;a2 = v2-((a1'*v2)/(a1'*a1))*a1;a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;A = [a1,a2,a3];计算得到的Rossler系统的LE为————0.0632310.092635-9.8924 Wolf文章中计算得到的Rossler系统的LE为————0.090-9.77需要注意的是——定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象。
chen系统李雅普诺夫指数的matlab源代码
chen系统李雅普诺夫指数的matlab源代码以下是计算 Chen 系统的李雅普诺夫指数的 Matlab 源代码。
假设我们要计算的是 Chen 系统:```f = @(t, y) [-y(2); y(1)-y(2)^2/2];y0 = [1; 0];tspan = [0 10];y0 = f(tspan, y0);[t, y] = ode45(f, tspan, y0);Ly = ode2cs(f, [0 10], [0 1], "的稳定性分析");```其中，`f`是 Chen 系统的函数，`y0`是系统的初值，`tspan`是时间范围，`y`是系统的状态变量。
`Ly`是李雅普诺夫指数，它表示系统的稳定性。
Matlab 代码的具体步骤如下:1. 定义 Chen 系统的函数`f`,并设置初值`y0`。
2. 使用`ode45`函数求解 Chen 系统的 ODE。
3. 使用`ode2cs`函数计算李雅普诺夫指数`Ly`。
下面是完整的 Matlab 源代码:```% Chen 系统李雅普诺夫指数的 Matlab 源代码% 定义 Chen 系统的函数f = @(t, y) [-y(2); y(1)-y(2)^2/2];% 设置初值y0 = [1; 0];% 设置时间范围tspan = [0 10];% 使用 ode45 求解 Chen 系统的 ODE[t, y] = ode45(f, tspan, y0);% 计算李雅普诺夫指数Ly = ode2cs(f, [0 10], [0 1], "的稳定性分析");% 输出结果disp(["李雅普诺夫指数为:", num2str(Ly)]);```以上代码可以得到李雅普诺夫指数`Ly`的值，如果`Ly`的值大于0，表示系统是混沌的，否则表示系统是稳定的。
rossler混沌系统lyapunov指数matlab实现代码[精彩]
![rossler混沌系统lyapunov指数matlab实现代码[精彩]](https://img.taocdn.com/s3/m/ba500cd32dc58bd63186bceb19e8b8f67c1cef1e.png)
Rossler 混沌系统Lyapunov指数Matlab实现代码Rossler 混沌系统Lyapunov指数Matlab实现代码分类：转载2010-04-11 10:49 混沌系统的基本特点就是系统对初始值的极端敏感性，两个相差无几的初值所产生的轨迹，随着时间的推移按指数方式分离，lyapunov指数就是定量的描述这一现象的量。
Lyapunov指数是衡量系统动力学特性的一个重要定量指标,它表征了系统在相空间中相邻轨道间收敛或发散的平均指数率。
对于系统是否存在动力学混沌, 可以从最大Lyapunov指数是否大于零非常直观的判断出来: 一个正的Lyapunov指数,意味着在系统相空间中,无论初始两条轨线的间距多么小,其差别都会随着时间的演化而成指数率的增加以致达到无法预测,这就是混沌现象。
Lyapunov指数的和表征了椭球体积的增长率或减小率,对Hamilton 系统,Lyapunov指数的和为零; 对耗散系统,Lyapunov指数的和为负。
如果耗散系统的吸引子是一个不动点,那么所有的Lyapunov指数通常是负的。
如果是一个简单的m维流形(m = 1或m = 2分别为一个曲线或一个面) ,那么,前m 个Lyapunov指数是零,其余的Lyapunov指数为负。
不管系统是不是耗散的,只要λ1 > 0就会出现混沌。
微分动力系统L yapunov指数的性质对于一维(单变量) 情形,吸引子只可能是不动点(稳定定态) 。
此时λ是负的。
对于二维情形, 吸引子或者是不动点或者是极限环。
对于不动点,任意方向的δxi , 都要收缩, 故这时两个Lyapunov指数都应该是负的, 即对于不动点, (λ1 ,λ 2 ) = ( - , - ) 。
至于极限环,如果取δxi 始终是垂直于环线的方向,它一定要收缩,此时λ < 0;当取δxi沿轨道切线方向,它既不增大也不缩小,可以想像,这时λ = 0。
事实上,所有不终止于定点而又有界的轨道(或吸引子) 都至少有一个Lyapunov指数等于零,它表示沿轨线的切线方向既无扩展又无收缩的趋势。
matlab练习程序（解西尔维斯特、李雅普诺夫方程）
matlab练习程序（解西尔维斯特、李雅普诺夫⽅程）西尔维斯特⽅程的形式：AX+XB=C李雅普诺夫⽅程的形式：AX+XA'=-C这两种⽅程都是已知矩阵A,B,C,求解X的⽅程。
对于这种⽅程有两种⽅法来求解，⼀种是朴素法，⼀种是Bartels-Stewart法。
以西尔维斯特⽅程为例，朴素法是将⽅程写为下列形式进⾏直接求解：其中圆圈中带个X的符号是kron积，vec()是将X或C转换为了n*1的列向量。
该⽅法将原来O(n^3)的问题变为了O(n^6)，如果矩阵⽐较⼤，估计速度会⽐较慢。
第⼆种⽅法为Bartels-Stewart法，下⾯以西尔维斯特⽅程为例介绍⼀下：⾸先我们对A和B‘进⾏shur分解：原⽅程可改写为：此时令：得到：此时R和S都是⼀个上三⾓矩阵，我们需要S作为下三⾓矩阵才能⽅便求解。
这⾥的S正好是我们是对B'的shur分解，由于shur分解的特性，shur(B)=VSV'，shur(B')=VS'V'，所以这⾥再对S求个转置即可。
得到类似下⾯的矩阵⽅程：展开得到：再依次求出y4,y3,y2,y1即可。
matlab代码如下：解西尔维斯特⽅程：clear all;close all;clc;A = rand(2,2);X = rand(2,2);B = rand(2,2);C = A*X+X*B;X%%系统函数X = sylvester(A,B,C)%%朴素法，⾃写kron积IA = [A zeros(2);zeros(2) A];BI = [B(1,1)*eye(2) B(2,1)*eye(2);B(1,2)*eye(2) B(2,2)*eye(2)];X = reshape(inv(IA+BI)*C(:),[2,2])%%朴素法，系统kron积X = reshape(inv(kron(eye(2),A)+kron(B',eye(2)))*C(:),[2,2])%%Bartels–Stewart法[U,R] = schur(A); %schur分解，R是上三⾓[V,S] = schur(B');S = S'; %S是下三⾓F = U'*C*V;%解R*Y + Y*S = F⽅程Y = zeros(2,2);Y(2,2) = F(2,2)/(R(2,2) + S(2,2));Y(2,1) = (F(2,1) - S(2,1)*Y(2,2)) / (S(1,1) + R(2,2));Y(1,2) = (F(1,2) - R(1,2)*Y(2,2)) / (R(1,1) + S(2,2));Y(1,1) = (F(1,1) - R(1,2)*Y(2,1) - S(2,1)*Y(1,2)) / (R(1,1) + S(1,1)); X = U*Y*V'解李雅普诺夫⽅程：clear all;close all;clc;A = rand(2);X = rand(2);C = -(A*X+X*A');X%%系统函数X = lyap(A,C)%%朴素法，⾃写kron积IA = [A zeros(2);zeros(2) A];BI = [A(1,1)*eye(2) A(1,2)*eye(2);A(2,1)*eye(2) A(2,2)*eye(2)];X = reshape(inv(IA+BI)*(-C(:)),[2,2])%%朴素法，系统kron积X = reshape(inv(kron(eye(2),A)+kron(A,eye(2)))*(-C(:)),[2,2]) %%Bartels–Stewart法[U,R] = schur(A); %schur分解，R是上三⾓S = R'; %S是下三⾓F = U'*(-C)*U;%解R*Y + Y*S = F⽅程Y = zeros(2,2);Y(2,2) = F(2,2)/(R(2,2) + S(2,2));Y(2,1) = (F(2,1) - S(2,1)*Y(2,2)) / (S(1,1) + R(2,2));Y(1,2) = (F(1,2) - R(1,2)*Y(2,2)) / (R(1,1) + S(2,2));Y(1,1) = (F(1,1) - R(1,2)*Y(2,1) - S(2,1)*Y(1,2)) / (R(1,1) + S(1,1)); X = U*Y*U'。
Lyapunov指数的计算方法
1【总结】Lyapunov 指数的计算方法非线性理论 近期为了把计算LE 的一些问题弄清楚，看了有7〜9本书！下面以吕金虎《混沌 时间序列分析及其应用》、马军海《复杂非线性系统的重构技术》为主线，把目前 已有的LE 计算方法做一个汇总！1.关于连续系统Lyapunov 指数的计算方法 连续系统LE 的计算方法主要有定义 方法、Jacobian 方法、QR 分解方法、奇异值分解方法，或者通过求解系统的微分 方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的 LE 求解 方法来计算得到。
关于连续系统 LE 的计算，主要以定义方法、Jacobian 方法做主 要介绍内容。
（1）定义法— 对H 堆连续动力系統z = 在—OBJ 孙 味“为中心.|拆（心0）||为丰笹啟存 «堆的球面*施著时间的演化，在t 时討该球而0P 变形为M 继的椭球厨・设滾椭域面的第/ 个坐标轴方向的半轴长対卩兀|，则诙系统第i 个指数対*此即连续系统Lyapunov 揩较飽定冥・弼计尊时・取|处（心0）［为岀W 为常数），以孔为球心・欧几里篇范敢为山的正衮 矢量集仙测，…叮为初始球.由非线性徴分方崔“尸㈤可以分别计算出点細 血创、 引他、r 引址经过时间t 后淺化的轨迹・役其终了点分别为珊、砒、f 仙 则令石f 陶一心■处严=甩-和，r 亦耳国=略一報#则可得新的矢重棄 叶禺巴…后畀｝・由于各牛妥臺在演化过程中舌焙向着是大的UapurOT IS 数方何靠掘，因此需要通过Schimdt IE 交化不断地讨新矢量逬行置换.SP Wolf to 文章中提出的GSR^法.表述如下：播着以他为球心，疤数対（I 的正奁矢臺料创巴叫叫…伽严;为祈球继续进行演 出 设演化至N步时，得到矢董慕冈㈤出巴…耳僅牛且N足够大,这可以得到Lyapunov 扌鐵的近似计算公式三实际计算时，取为1・定义法求解 Lyapunov 指数 JPG关于定义法求解的程序，和 matlab板块的 连续系统LE求解程序”差不多。
Matlab方程的求解
Matlab的Lyapunov、Sylvester和Riccati方程的Matlab求解一、连续Lyapunov方程连续Lyapunov方程可以表示为Lyapunov方程来源与微分方程稳定性理论，其中要求C为对称正定的n×n方阵，从而可以证明解X亦为n×n对称矩阵，这类方程直接求解比较困难，不过有了Matlab中lyap()函数，就简单多了。
>> A=[1 2 3;4 5 6;7 8 0]A =1 2 34 5 67 8 0>> C=-[10 5 4;5 6 7;4 7 9]C =-10 -5 -4-5 -6 -7-4 -7 -9>> X=lyap(A,C)X =-3.9444 3.8889 0.38893.8889 -2.7778 0.22220.3889 0.2222 -0.1111二、Lyapunov方程的解析解利用Kroncecker乘积的表示方法，可以将Lyapunov方程写为function x=lyap2(A,C)%Lyapunov方程的符号解法n=size(C,1);A0=kron(A,eye(n))+kron(eye(n),A);c=C(:);x0=-inv(A0)*c;x=reshape(x0,n,n)例子>>A=[1 2 3;4 5 6;7 8 0];>>C=-[10 5 4;5 6 7;4 7 9];>>x=lyap2(sym(A),sym(C))x =[ -71/18, 35/9, 7/18][ 35/9, -25/9, 2/9][ 7/18, 2/9, -1/9]三、离散Lyapunov方程离散Lyapunov方程的一般形式为Matlab中直接提供了dlyap()函数求解该方程：X=dlyap(A,Q)其实，如果A矩阵非奇异，则等式两边同时右乘得到就可以将其变换成连续的Sylvester方程而Sylvester方程是广义Lyapunov方程，故离散的Lyapunov方程还可以使用下面的方法求解B=-inv(A’)C=Q*in v(A’)X=lyap(A,B,C)下面总结下我们上面的讲到的知识点：X=lyap(A,C) 连续Lyapunov方程数值解法X=lyap2(A,C) 连续Lyapunov方程符号解法X=lyap(A,B,C) 广义Lyapunov方程，即Sylvester方程X=dlyap(A,Q)或者X=lyap(A,-inv(A’),Q*inv(A’))离散Lyapunov方程Sylvester方程Matlab求解Sylvester方程的一般形式为该方程又称为广义的Lyapunov方程，式中A是n×n方阵，B是m×m方阵，X 和C是n×m矩阵。
matlab编写的Lyapunov指数计算程序
% | y x -b |
%
% Then, the variational equation has a form:
%
% F = J*Y
% where Y is a square matrix with the same dimension as J.
% Corresponding m-file:
for l=1:n1 gsc(k)=gsc(k)+y(n1*j+l)*y(n1*k+l); end;
end;
for k=1:n1
for l=1:(j-1)
y(n1*j+k)=y(n1*j+k)-gsc(l)*y(n1*l+k);
end;
end;
znorm(j)=0.0;
for k=1:n1 znorm(j)=znorm(j)+y(n1*j+k)^2; end;
% stept - step on t-variable for Gram-Schmidt renormalization procedure.
% tend - finish value of time
% ystart - start point of trajectory of ODE system.
while calculation_progress == 1
[T,Y] = integrator(DS(1).method_int,@ode_lin,[t tt],y,options,P,n,neq,n_exp);
first_call = 0;
if calculation_progress == 99, break; end;
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
求解系统的Lyapunov指数谱程序
Lyapunov 指数是描述时序数据所生成的相空间中两个极其相近的初值所产生的轨道，随时间推移按指数方式分散或收敛的平均变化率。
任何一个系统，只要有一个Lyapunov 大于零，就认为该系统为混沌系统。
李雅普诺夫指数是指在相空间中相互靠近的两条轨线随着时间的推移，按指数分离或聚合的平均变化速率。
一 chen系统的Lyapunov指数谱
function dX = Chen2(t,X)
% Chen吸引子，用来计算Lyapunov指数
% dx/dt=a*(y-x)
% dy/dt=(c-a)*x+c*y-x*z
% dz/dt=x*y-b*z
global a; % 变量不放入参数表中
global b;
global c;
x=X(1); y=X(2); z=X(3);
% Y的三个列向量为相互正交的单位向量
Y = [X(4), X(7), X(10);
X(5), X(8), X(11);
X(6), X(9), X(12)];
% 输出向量的初始化
dX = zeros(12,1);
% Lorenz吸引子
dX(1) = a*(y-x);
dX(2) = (c-a)*x+c*y-x*z;
dX(3) = x*y-b*z;
% Lorenz吸引子的Jacobi矩阵
Jaco = [-a a 0;
c-a-z c -x;
y x -b];
dX(4:12) = Jaco*Y;
Z1=[];
Z2=[];
Z3=[];
global a;
global b;
global c;
b=3;c=28;
for a=linspace(32,40,100);
y=[1;1;1;1;0;0;0;1;0;0;0;1];
lp=0;
for k=1:200
[T,Y] = ode45('Chen2', 1, y);
y = Y(size(Y,1),:);
y0 = [y(4) y(7) y(10);
y(5) y(8) y(11);
y(6) y(9) y(12)];
y0=GS(y0);
mod(1)=norm(y0(:,1));
mod(2)=norm(y0(:,2));
mod(3)=norm(y0(:,3));
lp = lp+log(abs(mod));
y0(:,1)=y0(:,1)/mod(1);
y0(:,2)=y0(:,2)/mod(2);
y0(:,3)=y0(:,3)/mod(3);
y(4:12) = y0';
end
lp=lp/200;
Z1=[Z1 lp(1)];
Z2=[Z2 lp(2)];
Z3=[Z3 lp(3)];
end
a=linspace(32,40,100);
plot(a,Z1,'-',a,Z2,'-',a,Z3,'-');
title('Lyapunov exponents of Chen')
xlabel('b=3,c=28,parameter a'),ylabel('lyapunov exponents') grid on
以上是三个变量的Lyapunov指数谱，下面是最大的Lyapunov指数谱：
Z=[];
d0=1e-8;
for a=linspace(32,40,80)
lsum=0;
x=1;y=1;z=1;
x1=1;y1=1;z1=1+d0;
for i=1:100
[T1,Y1]=ode45('Chen',1,[x;y;z;a;3;28]);
[T2,Y2]=ode45('Chen',1,[x1;y1;z1;a;3;28]);
n1=length(Y1);n2=length(Y2);
x=Y1(n1,1);y=Y1(n1,2);z=Y1(n1,3);
x1=Y2(n2,1);y1=Y2(n2,2);z1=Y2(n2,3);
d1=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2);
x1=x+(d0/d1)*(x1-x);
y1=y+(d0/d1)*(y1-y);
z1=z+(d0/d1)*(z1-z);
if i>50
lsum=lsum+log(d1/d0);
end
end
Z=[Z lsum/(i-50)];
end
a=linspace(32,40,80);
plot(a,Z,'-');
title('Chen 系统最大lyapunov指数')
xlabel('parameter a'),ylabel('lyapunov exponents')
二模拟 Lorenz 系统最大lyapunov指数谱
function ly=jose_ly(b,k)
% the largest lyapunov exponent of josephson
% k 迭代步数，b 参数
% 方程如下:
% θ''+G*θ'+sinθ=I+A*sin(ωt)+αsin(βωt) % 变化：
% dx=y
% dy=-G*y-sin(x)+I+A*sin(w*t)+a*sin(b*w*t) %
% Example:
% ly=jose_ly(0,800)
%
% Author:LDYU
% Author's email: ustb03-
%
d0=1e-8;
ly=0;
lsum=0;
x=[0;2;b];
x1=[d0;2;b];
for t=1:k
[T1,Y1]=ode45('Josephon',[t-1,t],x);
[T2,Y2]=ode45('Josephon',[t-1,t],x1);
x=Y1(end,:);
x1=Y2(end,:);
d1=norm(x-x1);
x1=x+(d0/d1)*(x1-x);
lsum=lsum+log(d1/d0);
end
ly=lsum/k;。