Matlab画Lorenz系统的最大李雅普诺夫指数图
-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需要注意的是——定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象。
李雅普诺夫稳定性的Matlab实现
1 -1 -1
P -1
3
2
-1 2 5
对称矩阵的定号性(正定性)的判定(6/12)
Matlab程序m5-1如下。
P=[1 -1 -1; -1 3 2; -1 2 -5];
1 -1 -1
P -1 3
2
-1 2 5
result_state=posit_def(P); % 采用合同变换法判定矩阵定号性 switch result_state(1:5) % 运用开关语句,分类陈述矩阵正定否的判定结果
对称矩阵的定号性(正定性)的判定(9/12)
函数max()的主要调用格式为: [C,I]=max(A) D=max(A,B)
对第1种调用格式,若输入A为向量,输出的C为向量A的各 元素中最大值,输出I为该最大值在向量中的位置; 若A为矩阵,则C为矩阵A的各列的各元素中最大值,输 出I为这些最大值在各列的位置,这里输出C和I均为1 维数组。 如,执行语句 [C,I]=max([1 -2 3; -4 5 -6]); 后,C和I分别为[1 5 3]和[1 2 1]。
对称矩阵的定号性(正定性)的判定(10/12)
对第2种调用格式的输入A和B须为维数大小相同的矩阵 或向量,输出D为A和B两矩阵同样位置的元素的最大值组 成的矩阵。 如,执行语句 C=max([1 -2 3; -4 5 -6], [-1 2 -3; 4 -5 6]); 后,C为如下矩阵
1 2 3 4 5 6
Ch.5 李雅普诺夫稳定性 分析
目录
概述 5.1 李雅普诺夫稳定性的定义 5.2 李雅普诺夫稳定性的基本定理 5.3 线性系统的稳定性分析 5.4 非线性系统的稳定性分析 5.5 Matlab问题 本章小结
-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需要注意的是——定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象。
Lorenz模型论文稳定性论文：一个Lorenz模型的数值解法
Lorenz模型论文稳定性论文：一个Lorenz模型的数值解法摘要：lorenz方程是描述混沌现象的第一例有名的方程。
本文从气象学和解的稳定性两方面对lorenz模型进行概述，然后对此模型作了数值分析，基于稳定性讨论了程序设计并用matlab语言编写程序进行了求解，对数值解在图形上进行了描绘，最后对数值解的收敛性及稳定性进行了分析。
关键词：lorenz模型数值解收敛性稳定性第一章 lorenz模型概述本文要研究的lorenz方程形式如下：其中参数a=10,b=,c=28，初值条件为。
当时，原点是lorenz方程的唯一平衡点。
取李雅普诺夫函数，容易验证，因定正、定负，原点是渐进稳定的。
lorenz方程的所有轨迹均趋于原点。
当时，原点仍是lorenz方程的唯一平衡点。
但，，，出现叉式分支，原点不稳定。
而当时，原点不稳定。
且此时除原点外还出现两个异于原点的平衡点（，，），（,,），对称于轴。
对此两平衡点，考虑在平衡点处线性化lorenz方程，可求得特征方程为。
因，特征方程系数均大于0，实特征根必为负根。
知平衡点，渐近稳定的条件是，或，，其中。
当时，，特征方程有一对共轭纯虚根，出现分支；当时，特征方程除一负根外有一对共轭复特征根，其实部为正，对空间线性微分方程，这种空间平衡点为鞍焦点，空间轨迹投映于平面上为焦点和鞍点状。
固定a=10,b=进行讨论，易知此时=24.7368。
当时，lorenz方程的所有轨线趋于原点；当时，存在原点和平面上三个平衡点。
当时，平衡点是稳定的；当时，平衡点不稳定，属鞍焦点。
因为此时=28，所以平衡点不稳定，属鞍焦点。
取参数的不同值，我们可以通过数值解画出lorenz方程在相空间的轨迹图貌。
当时，由原点出发的两条轨迹各自分别趋于两平衡点，；在处，出现同宿轨；当时，出现由原点出发的两条轨线各自分别绕过一平衡点趋于另一平衡点，并在相空间中可能存在闭轨线或其他复杂轨线；当时，由于两平衡点，属鞍焦点，相空间中的轨线更为复杂。
洛伦兹系统李雅普诺夫指数的MATLAB源代码
%
% ****************************************************
% * 调用龙格库塔迭代公式进行迭代计算，算出下一时刻的值
% ****************************************************
lamda2 = lsum2/count1/h0;
lamda3 = lsum3/count1/h0;
% Lyapunov = [lamda1, lamda2, lamda3];
fprintf('LE1=%f,LE2=%f,LE3=%f\n',lamda1,lamda2,lamda3);
x3new = x3+dx3; y3new = y3+dy3; z3new = z3+dz3;
%
% ****************************************************
% * 求出前后两个时候的距离差值f
% ****************************************************
if j>= timesum1
lsum1 = lsum1+log(d1/d0);
lsum2 = lsum2+log(d2/d0);
lsum3 = lsum3+log(d3/d0);
count1 = count1+1;
%
f1x = x1new-x0new; f1y = y1new-y0new; f1z = z1new-z0new;
Lyapunov指数的计算方法
总结Lyapunov 指数的计算方法近期为了把计算LE的一些问题弄清楚,看了有7～9本书下面以吕金虎混沌时间序列分析及其应用、马军海复杂非线性系统的重构技术为主线,把目前已有的LE计算方法做一个汇总1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法,或者通过求解系统的微分方程,得到微分方程解的时间序列,然后利用时间序列即离散系统的LE求解方法来计算得到;关于连续系统LE的计算,主要以定义方法、Jacobian方法做主要介绍内容;1定义法关于定义法求解的程序,和matlab板块的“连续系统LE求解程序”差不多;以Rossler 系统为例Rossler系统微分方程定义程序function dX = Rossler_lyt,X% Rossler吸引子,用来计算Lyapunov指数% a=,b=,c=% dx/dt = -y-z,% dy/dt = x+ay,% dz/dt = b+zx-c,a = ;b = ;c = ;x=X1; y=X2; z=X3;% Y的三个列向量为相互正交的单位向量Y = X4, X7, X10;X5, X8, X11;X6, X9, X12;% 输出向量的初始化,必不可少dX = zeros12,1;% Rossler吸引子dX1 = -y-z;dX2 = x+ay;dX3 = b+zx-c;% Rossler吸引子的Jacobi矩阵Jaco = 0 -1 -1;1 a 0;z 0 x-c;dX4:12 = JacoY;求解LE代码：% 计算Rossler吸引子的Lyapunov指数clear;yinit = 1,1,1;orthmatrix = 1 0 0;0 1 0;0 0 1;a = ;b = ;c = ;y = zeros12,1;% 初始化输入y1:3 = yinit;y4:12 = orthmatrix;tstart = 0; % 时间初始值tstep = 1e-3; % 时间步长wholetimes = 1e5; % 总的循环次数steps = 10; % 每次演化的步数iteratetimes = wholetimes/steps; % 演化的次数mod = zeros3,1;lp = zeros3,1;% 初始化三个Lyapunov指数Lyapunov1 = zerositeratetimes,1;Lyapunov2 = zerositeratetimes,1;Lyapunov3 = zerositeratetimes,1;for i=1:iteratetimestspan = tstart:tstep:tstart + tstepsteps;T,Y = ode45'Rossler_ly', tspan, y;% 取积分得到的最后一个时刻的值y = YsizeY,1,:;% 重新定义起始时刻tstart = tstart + tstepsteps;y0 = y4 y7 y10;y5 y8 y11;y6 y9 y12;%正交化y0 = ThreeGSy0;% 取三个向量的模mod1 = sqrty0:,1'y0:,1;mod2 = sqrty0:,2'y0:,2;mod3 = sqrty0:,3'y0:,3;y0:,1 = y0:,1/mod1;y0:,2 = y0:,2/mod2;y0:,3 = y0:,3/mod3;lp = lp+logabsmod;%三个Lyapunov指数Lyapunov1i = lp1/tstart;Lyapunov2i = lp2/tstart;Lyapunov3i = lp3/tstart;y4:12 = y0';end% 作Lyapunov指数谱图i = 1:iteratetimes;ploti,Lyapunov1,i,Lyapunov2,i,Lyapunov3程序中用到的ThreeGS程序如下：%G-S正交化function A = ThreeGSV % V 为33向量v1 = V:,1;v2 = V:,2;v3 = V:,3;a1 = zeros3,1;a2 = zeros3,1;a3 = zeros3,1;a1 = v1;a2 = v2-a1'v2/a1'a1a1;a3 = v3-a1'v3/a1'a1a1-a2'v3/a2'a2a2;A = a1,a2,a3;计算得到的Rossler系统的LE为————Wolf文章中计算得到的Rossler系统的LE为———— 0需要注意的是——定义法求解的精度有限,对有些系统的计算往往出现计果和理论值有偏差的现象;正交化程序可以根据上面的扩展到NN向量,这里就不加以说明了,对matlab用户来说应该还是比较简单的2Jacobian方法通过资料检索,发现论坛中用的较多的LET工具箱的算法原理就是Jacobian方法;基本原理就是首先求解出连续系统微分方程的近似解,然后对系统的Jacobian矩阵进行QR分解,计算Jacobian矩阵特征值的乘积,最后计算出LE和分数维;经过计算也证明了这种方法精度较高,对目前常见的混沌系统,如Lorenz、Henon、Duffing等的Lyapunov指数的计算精度都很高,而且程序编写有一定的规范,个人很推荐使用;虽然我自己要做的系统并不适用LET工具箱可以在网络上找到,这里就不列出了关于LET工具箱如果有问题,欢迎加入本帖讨论Jacobian法求解Lyapunov指数.JPG对离散动力系统,或者说是非线性时间序列,往往不需要计算出所有的Lyapunov指数,通常只需计算出其最大的Lyapunov指数即可;“1983年,格里波基证明了只要最大Lyapunov指数大于零,就可以肯定混沌的存在”;目前常用的计算混沌序列最大Lyapunov指数的方法主要有一下几种：1由定义法延伸的Nicolis方法2Jacobian方法3Wolf方法4P－范数方法5小数据量方法其中以Wolf方法和小数据量方法应用最为广泛,也最为普遍;下面对Nicolis方法、Wolf方法以及小数据量方法作一一介绍;1Nicolis方法这种方法和连续系统的定义方法类似,而且目前应用很有限制,因此只对其理论进行介绍,编程应用方面就省略了Nicolis方法求最大2Wolf方法Wolf方法求最大Wolf方法的Matlab程序如下：function lambda_1=lyapunov_wolfdata,N,m,tau,P% 该函数用来计算时间序列的最大Lyapunov 指数--Wolf 方法% m: 嵌入维数一般选大于等于10% tau:时间延迟一般选与周期相当,如我选2000% data:时间序列可以选1000；% N:时间序列长度满足公式：M=N-m-1tau=24000-10-11000=5000% P:时间序列的平均周期,选择演化相点距当前点的位置差,即若当前相点为I,则演化相点只能在|I－J|>P的相点中搜寻P=周期=600% lambda_1: 返回最大lyapunov指数值min_point=1 ; %&&要求最少搜索到的点数MAX_CISHU=5 ; %&&最大增加搜索范围次数%FLYINGHAWK% 求最大、最小和平均相点距离max_d = 0; %最大相点距离min_d = +100; %最小相点距离avg_dd = 0;Y=reconstitutiondata,N,m,tau; %相空间重构可将此段程序加到整个程序中,在时间循环内,可以保存时间序列的地方;见完整程序;M=N-m-1tau; %重构相空间中相点的个数for i = 1 : M-1for j = i+1 : Md = 0;for k = 1 : md = d + Yk,i-Yk,jYk,i-Yk,j;endd = sqrtd;if max_d < dmax_d = d;endif min_d > dmin_d = d;endavg_dd = avg_dd + d;endendavg_d = 2avg_dd/MM-1; %平均相点距离dlt_eps = avg_d - min_d ; %若在min_eps～max_eps中找不到演化相点时,对max_eps的放宽幅度min_eps = min_d + dlt_eps / 2 ; %演化相点与当前相点距离的最小限max_eps = min_d + 2 dlt_eps ; %&&演化相点与当前相点距离的最大限% 从P+1～M-1个相点中找与第一个相点最近的相点位置Loc_DK及其最短距离DKDK = +100; %第i个相点到其最近距离点的距离Loc_DK = 2; %第i个相点对应的最近距离点的下标for i = P+1:M-1 %限制短暂分离,从点P+1开始搜索d = 0;for k = 1 : md = d + Yk,i-Yk,1Yk,i-Yk,1;endd = sqrtd;if d < DK & d > min_epsDK = d;Loc_DK = i;endend% 以下计算各相点对应的李氏数保存到lmd数组中% i 为相点序号,从1到M-1,也是i-1点的演化点；Loc_DK为相点i-1对应最短距离的相点位置,DK为其对应的最短距离% Loc_DK+1为Loc_DK的演化点,DK1为i点到Loc_DK+1点的距离,称为演化距离% 前i个log2DK1/DK的累计和用于求i点的lambda值sum_lmd = 0 ; % 存放前i个log2DK1/DK的累计和for i = 2 : M-1 % 计算演化距离DK1 = 0;for k = 1 : mDK1 = DK1 + Yk,i-Yk,Loc_DK+1Yk,i-Yk,Loc_DK+1;endDK1 = sqrtDK1;old_Loc_DK = Loc_DK ; % 保存原最近位置相点old_DK=DK;% 计算前i个log2DK1/DK的累计和以及保存i点的李氏指数if DK1 ~= 0& DK ~= 0sum_lmd = sum_lmd + logDK1/DK /log2;endlmdi-1 = sum_lmd/i-1; 此处可以保存不同相点i对应的李氏指数,见完整程序;; % 以下寻找i点的最短距离:要求距离在指定距离范围内尽量短,与DK1的角度最小 point_num = 0; % &&在指定距离范围内找到的候选相点的个数cos_sita = 0 ; %&&夹角余弦的比较初值——要求一定是锐角zjfwcs=0 ; %&&增加范围次数while point_num == 0% 搜索相点for j = 1 : M-1if absj-i <=P-1 %&&候选点距当前点太近,跳过continue;end%计算候选点与当前点的距离dnew = 0;for k = 1 : mdnew = dnew + Yk,i-Yk,jYk,i-Yk,j;enddnew = sqrtdnew;if dnew < min_eps| dnew > max_eps %&&不在距离范围,跳过continue;end%计算夹角余弦及比较DOT = 0;for k = 1 : mDOT = DOT+Yk,i-Yk,jYk,i-Yk,old_Loc_DK+1;endCTH = DOT/dnewDK1;if acosCTH > 3./4 %&&不是小于45度的角,跳过continue;endif CTH > cos_sita %&&新夹角小于过去已找到的相点的夹角,保留cos_sita = CTH;Loc_DK = j;DK = dnew;endpoint_num = point_num +1;endif point_num <= min_pointmax_eps = max_eps + dlt_eps;zjfwcs =zjfwcs +1;if zjfwcs > MAX_CISHU %&&超过最大放宽次数,改找最近的点DK = +100;for ii = 1 : M-1if absi-ii <= P-1 %&&候选点距当前点太近,跳过continue;endd = 0;for k = 1 : md = d + Yk,i-Yk,iiYk,i-Yk,ii;endd = sqrtd;if d < DK & d > min_epsDK = d;Loc_DK = ii;endendbreak;endpoint_num = 0 ; %&&扩大距离范围后重新搜索cos_sita = 0;endendend%取平均得到最大李雅普诺夫指数此处只有一个值,若为正说明体系是混沌的,若为负则说明体系是周期性的确定性运动lambda_1=sumlmd/lengthlmd;程序中用到的reconstitution函数如下：此段程序可直接放在时间循环内部,即可以保存时间序列的地方;见完整程序范例;function X=reconstitutiondata,N,m,tau%该函数用来重构相空间% m 为嵌入空间维数% tau 为时间延迟% data 为输入时间序列% N 为时间序列长度% X 为输出,是mn维矩阵M=N-m-1tau; %相空间中点的个数for j=1:M %相空间重构for i=1:mXi,j=datai-1tau+j;endend这里声明一下,这些程序并非我自己编写的,均是转载,其使用我已经验证过,绝对可以运行3小数据量方法说小数据量方法是目前最实用、应用最广泛的方法应该不为过吧,呵呵小数据量方法求最大Lyapunov指数.JPG上面两种方法,即Wolf方法和小数据量方法,在计算LE之前,都要求对时间序列进行重构相空间,重构相空间的优良对于最大LE的计算精度影响非常大因此重构相空间的几个参数的确定就非常重要;1时间延迟主要推荐两种方法——自相关函数法、C－C方法自相关函数法——对一个混沌时间序列,可以先写出其自相关函数,然后作出自相关函数关于时间t的函数图像;根据数值试验结果,当自相关函数下降到初始值的1－1/e时,所得的时间t即为重构相空间的时间延迟;C－C方法——可以同时计算出时间延迟和时间窗口,个人推荐使用这种方法2平均周期平均周期的计算可以采用FFT方法;在matlab帮助中有一个太阳黑子的例子,现摘录如下：load %装载数据文件year = sunspot:,1; %读取年份信息wolfer = sunspot:,2; %读取黑子活动数据plotyear,wolfer %绘制原始数据图title'Sunspot Data'Y = fftwolfer; %快速FFT变换N = lengthY; %FFT变换后数据长度Y1 = ; %去掉Y的第一个数据,它是所有数据的和power = absY1:N/2.^2; %求功率谱nyquist = 1/2;freq = 1:N/2/N/2nyquist; %求频率plotfreq,power, grid on %绘制功率谱图xlabel'cycles/year'title'Periodogram'period = 1./freq; %年份周期plotperiod,power, axis0 40 0 2e7, grid on ％绘制年份－功率谱曲线ylabel'Power'xlabel'PeriodYears/Cycle'mp,index = maxpower; %求最高谱线所对应的年份下标periodindex %由下标求出平均周期3嵌入维数目前嵌入维数的主要计算方法是采用Grassberger和Procaccia提出的G－P算法计算出序列的关联维数d,然后利用嵌入维数m>=2d+1,选取合适的嵌入维数;G—P算法程序如下：function ln_r,ln_C=G_Pdata,N,tau,min_m,max_m,ss% the function is used to calculate correlation dimention with G-P algorithm% 计算关联维数的G－P算法% data:the time series 时间序列% N: the length of the time series 时间序列长度% tau: the time delay 时间延迟% min_m:the least embedded dimention m 最小的嵌入维数% max_m:the largest embedded dimention m 最大的嵌入维数% ss:the stepsize of r r的步长%skyhawkfor m=min_m:max_mY=reconstitutiondata,N,m,tau;%reconstitute state spaceM=N-m-1tau;%the number of points in state spacefor i=1:M-1for j=i+1:Mdi,j=maxabsY:,i-Y:,j;%calculate the distance of each twoend %points in state space 计算状态空间中每两点之间的距离endmax_d=maxmaxd;%the max distance of all points 得到所有点之间的最大距离d1,1=max_d;min_d=minmind;%the min distance of all points 得到所有点间的最短距离 delt=max_d-min_d/ss;%the stepsize of r 得到r的步长for k=1:ssr=min_d+kdelt;Ck=correlation_integralY,M,r;%calculate the correlation integralln_Cm,k=logCk;%lnCrln_rm,k=logr;%lnrfprintf'%d/%d/%d/%d\n',k,ss,m,max_m;endplotln_rm,:,ln_Cm,:;hold on;endfid=fopen'','w';fprintffid,'% %\n',ln_r;fclosefid;fid = fopen'','w';fprintffid,'% %\n',ln_C;fclosefid;程序中的correlation_integral函数如下：function C_I=correlation_integralX,M,r%the function is used to calculate correlation integral%C_I:the value of the correlation integral%X:the reconstituted state space,M is a mM matrix%m:the embedding demention%M:M is the number of embedded points in m-dimensional sapce%r:the radius of the Heaviside function,sigma/2<r<2sigma%calculate the sum of all the values of Heaviside%skyhawksum_H=0;for i=1:M% fprintf'%d/%d\n',i,M;for j=i+1:Md=normX:,i-X:,j,inf;%calculat the distances of each two points in matris M with sup-normsita=heavisider,d;%calculate the value of the heaviside functionsum_H=sum_H+sita;endendC_I=2sum_H/MM-1;%the value of correlation integral发布于2007-08-03 20:35:01感谢octopussheng的总结这是现有的方法的一个总结,不知道你对这些方法有些什么体会,或者说他们的局限,现在还有作这方面的改进的吗先总结总结,具体的比较会逐步贴出来的,而且这也需要大家的一起努力啊光是我在这里说,也只是我自己的一点体会,希望用过这些方法的一起来参与这个总结工作以上的各种方法在实际应用的时候要根据具体情况来选择;一般地,如果已知系统方程当然系统不能太过复杂时,则计算Lyapunov指数采用定义法、Jacobian方法要精确、简单些而如果系统方程比较复杂如超维系统、或者为一时间序列时,则推荐采样Wolf 方法、小数据量方法;Wolf方法的特点是时间序列无噪声,空间中小向量的演变高度非线性,而Jacobian 方法则是噪声大,空间中小向量的演变接近线性;小数据量方法的优点在于：1对小数据组的计算可靠；2计算量较小,比wolf方法快很多；3编程、操作较为容易;而关于时间延迟、嵌入维数、平均周期的确定,还是推荐C－C方法和G－P算法,结果更为可靠一些你这已经做得很好了;这个具体的应用我也是做得很少的,还有请注意一下matlab版中的一个帖子, 早就注意了,我这里面一些代码就是从那边套用的试过了吗都没有问题吧计算起来,原来论坛里面一些程序不能运行,套用,那你就要注意,在这里重写就要有我们的特色;不能够是简单的重请详细读帖子,上面有说明,当然不是简单的重复啦,我这里是理论结合实践,呵呵程序都测试能算的,不过对不同的系统,精度不敢保证的太高,呵呵,只要能运行就可以参照自己来编写精度可以根据自己需要来选择呵呵,确实,保证能够使用希望这样总结一下对大家有所帮助吧非常有价值;谢谢。
求最大李雅普诺夫指数的matlab程序
程序一function dx=Lorenz(t,x);dx(1,1)=10*(x(2)-x(1));dx(2,1)=x(1)*(30-x(3))-x(2);dx(3,1)=x(1)*x(2)-8/3*x(3);dx(4,1)=0;dx(5,1)=0;dx(6,1)=0;function lambda_1=lyapunov_wolf1(data,N,m,tau,P)% 该函数用来计算时间序列的最大Lyapunov 指数--Wolf 方法% m: 嵌入维数% tau:时间延迟% data:时间序列% N:时间序列长度% P:时间序列的平均周期,选择演化相点距当前点的位置差，即若当前相点为I，则演化相点只能在|I－J|>P的相点中搜寻% lambda_1:返回最大lyapunov指数值%**************************************************************************% ode计算整数阶系统的时间序列%******************************************************************delt_t1 = 0.001;t1 = 0:delt_t1:60;[tt1,y1]=ode45(@lorenz,t1,[-1,0,1]);xx1 = y1(:,1)';x1 = spline(tt1, xx1, t1);data= x1(20000:10:60000);%采样N=length(data);m=3;tau=11;%*****************************************************% FFT计算平均周期%**********************************************************x=data;xPower=abs(fft(x)).^2;NN=length(xPower);xPower(1)=[];%去除直流分量NN=floor(NN/2);xPower=xPower(1:NN);freq=(1:NN)/NN*0.5;[mP,index]=max(xPower);P=index;%*************************************************************min_point=1 ; %&&要求最少搜索到的点数MAX_CISHU=5 ; %&&最大增加搜索范围次数%FL YINGHAWK% 求最大、最小和平均相点距离max_d = 0; %最大相点距离min_d = 1.0e+100; %最小相点距离avg_dd = 0;Y=reconstitution(data,N,m,tau); %相空间重构M=N-(m-1)*tau; %重构相空间中相点的个数for i = 1 : (M-1)for j = i+1 : Md = 0;for k = 1 : md = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));endd = sqrt(d);if max_d < dmax_d = d;endif min_d > dmin_d = d;endavg_dd = avg_dd + d;endendavg_d = 2*avg_dd/(M*(M-1)); %平均相点距离dlt_eps = (avg_d - min_d) * 0.02 ; %若在min_eps～max_eps中找不到演化相点时，对max_eps的放宽幅度min_eps = min_d + dlt_eps / 2 ; %演化相点与当前相点距离的最小限max_eps = min_d + 2 * dlt_eps ; %&&演化相点与当前相点距离的最大限% 从P+1～M-1个相点中找与第一个相点最近的相点位置(Loc_DK)及其最短距离DK DK = 1.0e+100; %第i个相点到其最近距离点的距离Loc_DK = 2; %第i个相点对应的最近距离点的下标for i = (P+1):(M-1) %限制短暂分离，从点P+1开始搜索d = 0;for k = 1 : md = d + (Y(k,i)-Y(k,1))*(Y(k,i)-Y(k,1));endd = sqrt(d);if (d < DK) & (d > min_eps)DK = d;Loc_DK = i;endend% 以下计算各相点对应的李氏数保存到lmd()数组中% i 为相点序号，从1到(M-1)，也是i-1点的演化点；Loc_DK为相点i-1对应最短距离的相点位置，DK为其对应的最短距离% Loc_DK+1为Loc_DK的演化点，DK1为i点到Loc_DK+1点的距离，称为演化距离% 前i个log2（DK1/DK）的累计和用于求i点的lambda值sum_lmd = 0 ; % 存放前i个log2（DK1/DK）的累计和for i = 2 : (M-1) % 计算演化距离DK1 = 0;for k = 1 : mDK1 = DK1 + (Y(k,i)-Y(k,Loc_DK+1))*(Y(k,i)-Y(k,Loc_DK+1));endDK1 = sqrt(DK1);old_Loc_DK = Loc_DK ; % 保存原最近位置相点old_DK=DK;% 计算前i个log2（DK1/DK）的累计和以及保存i点的李氏指数if (DK1 ~= 0)&( DK ~= 0)sum_lmd = sum_lmd + log(DK1/DK) /log(2);endlmd(i-1) = sum_lmd/(i-1);% 以下寻找i点的最短距离：要求距离在指定距离范围内尽量短，与DK1的角度最小 point_num = 0 ; % &&在指定距离范围内找到的候选相点的个数cos_sita = 0 ; %&&夹角余弦的比较初值——要求一定是锐角zjfwcs=0 ;%&&增加范围次数while (point_num == 0)% * 搜索相点for j = 1 : (M-1)if abs(j-i) <=(P-1) %&&候选点距当前点太近，跳过！continue;end%*计算候选点与当前点的距离dnew = 0;for k = 1 : mdnew = dnew + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));enddnew = sqrt(dnew);if (dnew < min_eps)|( dnew > max_eps ) %&&不在距离范围，跳过！continue;end%*计算夹角余弦及比较DOT = 0;for k = 1 : mDOT = DOT+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,old_Loc_DK+1));endCTH = DOT/(dnew*DK1);if acos(CTH) > (3.14151926/4) %&&不是小于45度的角，跳过！continue;endif CTH > cos_sita %&&新夹角小于过去已找到的相点的夹角，保留 cos_sita = CTH;Loc_DK = j;DK = dnew;endpoint_num = point_num +1;endif point_num <= min_pointmax_eps = max_eps + dlt_eps;zjfwcs =zjfwcs +1;if zjfwcs > MAX_CISHU %&&超过最大放宽次数，改找最近的点DK = 1.0e+100;for ii = 1 : (M-1)if abs(i-ii) <= (P-1) %&&候选点距当前点太近，跳过！continue;endd = 0;for k = 1 : md = d + (Y(k,i)-Y(k,ii))*(Y(k,i)-Y(k,ii));endd = sqrt(d);if (d < DK) & (d > min_eps)DK = d;Loc_DK = ii;endendbreak;endpoint_num = 0 ; %&&扩大距离范围后重新搜索cos_sita = 0;endendend%取平均得到最大李雅普诺夫指数lambda_1=sum(lmd)/length(lmd)function X=reconstitution(data,N,m,tau)%该函数用来重构相空间% m为嵌入空间维数% tau为时间延迟% data为输入时间序列% N为时间序列长度% X为输出,是m*n维矩阵M=N-(m-1)*tau;%相空间中点的个数for j=1:M %相空间重构for i=1:mX(i,j)=data((i-1)*tau+j);endend以上是计算最大李氏指数的程序，可以运行。
最大李雅谱诺夫数数值计算
最大李雅谱诺夫数数值计算————————————————————————————————作者：————————————————————————————————日期：最大李雅谱诺夫指数的数值计算------------------------------------------------------------------------最大李雅谱诺夫指数的数值计算J.C.SprottDepartment of Physics, University of Wisconsin, Madison, WI 53706, USAOctober 15, 1997(Revised August 31, 2004)袁来翻译/上海大学/05/03/22原文出处:/chaos/lyapexp.htm通常通过计算最大的李雅谱诺夫指数来检验系统是否混沌，一个正的最大李雅谱诺夫指数标致了该系统是否混沌，当你已经获得已知产生混沌的方程，这是相对容易做的。
当你仅仅只有试验数据记录时，这种计算的困难程度是几乎不可能的，这里我们不考虑这种情形。
计算的一般思路是跟踪两个较近的轨道，计算它们分离的平均对数率，无论它们分离得多远，其中之一的轨道将最终回到另一条轨道的附近。
在每次迭代过程中，一程序将实现这个工作，完整的程序过程如下：1：在吸引域内，开始于任何初值条件。
最好开始于吸引子上的一个已知点，这种情形，第2步可以省略掉。
2：迭代直至轨道落到吸引子上。
这需要自己的判断或通过研究得到的有关该系统的一些知识，对绝大多数系统，只要迭代几百次，我们就认为足够了。
通常情况下，除非你选的初值点刚好靠近分岔点，要不在几百次后它不会远离吸引子的。
3：选择（几乎任意的）靠近的点（设两点的距离为d0）适合的d0的选择是至少要比计算机中的浮点数的1000倍大，如：在（8字节）的双精度的（对这些计算的最低推荐），变量有52－比特的尾数，故，精度是：2-52=2.22x10-16，因此，我们只要有d0=10-12 就足够了。