时滞系统的最大lyapunov指数估算？
第二，如果不一样，那一般应用什么方法去求时滞系统的最大lyapunov指数？
两者没有关系，你给出的是运动微分方程，一般的处理思路是现对微分方程进行求解，求的方程的时间序列，然后再用你所说的wolf法等方法对相关序列进行分析求解lyapunov指数 你有什么好的带有程序的参考资料吗？ 我把写的程序发上来，但是有问题，请问怎么修改？function dx=loren(t,x,z)
%global r k4 b d1 a A B u mu1 C d2 tau
xlag=z(:,1);
s=zeros(4,1);
%% 模型形式
s(1)=50+1000.*x(2).*x(4)-0.000098.*x(3).*x(1)-0.011.*x(1);
s(2)=0.000098.*x(3).*x(1)-0.1.*x(2)-1200.*x(2).*x(4);
s(3)=200.*x(2)-0.67.*x(3);
s(4)=exp(-0.5).*0.05.*xlag(2).*xlag(4)-0.3.*x(4);
function lambda_1=lyapunov_wolf(data,N,m,tau,p)
%该函数用来计算时间序列的最大Lyapunov 指数--Wolf 方法
%m: 嵌入维数
%tau:时间延迟
%data:时间序列
%N:时间序列长度
%lambda_1:返回最大lyapunov指数值
%**************************************************************************
%ode计算整数阶系统的时间序列
%******************************************************************
Delays=;
delt_t1 = 0.001;
t1 = 0:delt_t1:60;
y1=dde23('loren',Delays,,t1);
xx1 = y1(:,1)';
x1 = spline(y1.x, xx1, t1);
data= x1(20000:10:60000); %采样
N=length(data);
m=4;
tau=3;
%*****************************************************
% 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;
=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; %第i个相点到其最近距离点的距离
Loc_DK = 2; %第i个相点对应的最近距离点的下标
for i = (P+1):(M-1) %限制短暂分离，从点P+1开始搜索
d = 0;
for k = 1 : m
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 ; % 存放前i个log2（DK1/DK）的累计和
for i = 2 : (M-1) % 计算演化距离
DK1 = 0;
for k = 1 : m
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) %&&不是小于45度的角，跳过！
continue;
end
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);
plot(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暂时运行结果图：
结果还是我自己搞出来啦。 1. wolf计算的最大lyapunov指数值，请问下图的横坐标是表示什么？
2. 下图t1的初值与末值应怎样取最好？这个是否求数值解的时间的初值与末值有关？
3. 为什么下图这个只取y1(:,1)'，而不取y1(:,2)'，y1(:,3)'？
可以分享一下吗，
页:
[1]