第1页 / 共6页
第2页 / 共6页
第3页 / 共6页
第4页 / 共6页
第5页 / 共6页
第6页 / 共6页
求最大李雅普诺夫指数的matlab程序.doc
程序一
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 ; %&&最大增加搜索范围次数
%FLYINGHAWK
% 求最大、最小和平均相点距离
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 : M
d = 0;
for k = 1 : m
%最大相点距离
%最小相点距离
%相空间重构
%重构相空间中相点的个数
d = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));
end
d = sqrt(d);
if max_d < d
max_d = d;
end
if min_d > d
min_d = d;
end
avg_dd = avg_dd + d;
end
end
avg_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;
Loc_DK = 2;
for i = (P+1):(M-1)
d = 0;
for k = 1 : m
%第 i 个相点到其最近距离点的距离
%第 i 个相点对应的最近距离点的下标
%限制短暂分离，从点 P+1 开始搜索
d = d + (Y(k,i)-Y(k,1))*(Y(k,i)-Y(k,1));
end
d = sqrt(d);
if (d < DK) & (d > min_eps)
DK = d;
Loc_DK = i;
end
end
% 以下计算各相点对应的李氏数保存到 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 ;
for i = 2 : (M-1)
DK1 = 0;
for k = 1 : m
% 存放前 i 个 log2（DK1/DK）的累计和
% 计算演化距离
DK1 = DK1 + (Y(k,i)-Y(k,Loc_DK+1))*(Y(k,i)-Y(k,Loc_DK+1));
end
DK1 = 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);
end
lmd(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 : m
dnew = dnew + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));
end
dnew = sqrt(dnew);
if (dnew < min_eps)|( dnew > max_eps ) %&&不在距离范围，跳过！
continue;
end
%*计算夹角余弦及比较
DOT = 0;
for k = 1 : m
DOT = DOT+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,old_Loc_DK+1));
end
CTH = DOT/(dnew*DK1);
if acos(CTH) > (3.14151926/4)
continue;
end
%&&不是小于 45 度的角，跳过！
if CTH > cos_sita %&&新夹角小于过去已找到的相点的夹角，保留
cos_sita = CTH;
Loc_DK = j;
DK = dnew;
end
point_num = point_num +1;
end
if point_num <= min_point
max_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;
end
d = 0;
for k = 1 : m
%&&候选点距当前点太近，跳过！
d = d + (Y(k,i)-Y(k,ii))*(Y(k,i)-Y(k,ii));
end
d = sqrt(d);
if (d < DK) & (d > min_eps)
DK = d;
Loc_DK = ii;
end
end
break;
end
point_num = 0
cos_sita = 0;
end
; %&&扩大距离范围后重新搜索
end
end
%取平均得到最大李雅普诺夫指数
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:m
X(i,j)=data((i-1)*tau+j);
end
end
以上是计算最大李氏指数的程序，可以运行。问题是里面的时间延迟 tau 和嵌入维数 m 该如
何取值，程序里我是随便取的 3 和 11
*******************************************************************************
程序二
现我做了一个基于 RHR 算法的李氏指数计算方法，收敛速度很快，精确度也还可
以，现传上：
Lya1=[];Lya2=[];Lya3=[];
V=eye(3);
S=V;b1=0;
a=0.4;c=0.2;gama=3.5;
b=4.0;
h=0.01;
x(1)=0.1;y(1)=0;z(1)=0;n=0;
while z<=200
n=n+1;
k1=h*y(n);
m1=h*(-sin(x(n))-a*y(n)+b*cos(gama*z(n)).*sin(x(n))+c);
k2=h*(y(n)+m1/2);
m2=h*(-sin(x(n)+k1/2)-a*(y(n)+m1/2)+b*cos(gama*(z(n)+h/2)).*sin(x(n)+k1/2)+c);
k3=h*(y(n)+m2/2);
m3=h*(-sin(x(n)+k2/2)-a*(y(n)+m2/2)+b*cos(gama*(z(n)+h/2)).*sin(x(n)+k2/2)+c);
k4=h*(y(n)+m3);
m4=h*(-sin(x(n)+k3)-a*(y(n)+m3)+b*cos(gama*(z(n)+h)).*sin(x(n)+k3)+c);
x(n+1)=x(n)+(k1+2*k2+2*k3+k4)/6;
y(n+1)=y(n)+(m1+2*m2+2*m3+m4)/6;
z(n+1)=n*h;
J = [0 1 0;
b*cos(gama*z(n+1))*cos(x(n+1))-cos(x(n+1)) -a -b*gama*sin(gama*z(n+1))*sin(x(n+1));
0 0 0];
J=eye(3)+h*J;
B=J*V*S;
[V,S,U]=svd(B);
a_max=max(diag(S));
S=(1/a_max)*S;
b1=b1+log(a_max);
Lyapunov1=(log(diag(S))+b1)/(n*h);
Lya1=[Lya1,Lyapunov1(1,
Lya2=[Lya2,Lyapunov1(2,
Lya3=[Lya3,Lyapunov1(3,
];
];
];
end
Lyapunov1
n=1:20001;
plot(n,Lya1,'k',n,Lya2,'k',n,Lya3,'k')
%grid on
axis([0,30001,-0.8,0.5])
title('Lyapunov exponents of Warship')
xlabel('n'),ylabel('LCE')
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
