how to find largest Lyapunov exponent for Lorenz in Rosenstein's Algorithm?

7 次查看（过去 30 天）
Syed Zubair shah
Syed Zubair shah 2022-5-5
回答： SOUMNATH PAUL 2024-1-17
function d = lyarosenstein(x,f,m,tao,meanperiod,maxiter)
% d:divergence of nearest trajectoires
% x:signal
% tao:time delay
% m:embedding dimension
N=length(x);
% [T, X] = feval(@ode45, @rossler, [1 100 + .01], x);
M=N-(m-1)*tao;
Y=psr_deneme(x,m,tao);
for i=1:M
x0=ones(M,1)*Y(i,:);
distance=sqrt(sum((Y-x0).^2,2));
for j=1:M
if abs(j-i)<=meanperiod
distance(j)=1e10;
end
end
[neardis(i) nearpos(i)]=min(distance);
end
for k=1:maxiter
maxind=M-k;
evolve=0;
pnt=0;
for j=1:M
if j<=maxind && nearpos(j)<=maxind
dist_k=sqrt(sum((Y(j+k,:)-Y(nearpos(j)+k,:)).^2,2));
if dist_k~=0
evolve=evolve+log(dist_k);
pnt=pnt+1;
end
end
end
if pnt > 0
d(k)=evolve/pnt;
else
d(k)=0;
end
end
figure
plot(log(d))
% LLE Calculation
fs=2000;%sampling frequency
tlinear=15:78;
F = polyfit(tlinear,d(tlinear),1);
lle = F(1)*fs;
function Y=psr_deneme(x,m,tao,npoint)
%Phase space reconstruction
%x : time series
%m : embedding dimension
%tao : time delay
%npoint : total number of reconstructed vectors
%Y : M x m matrix
% author:"Merve Kizilkaya"
N=length(x);
if nargin == 4
M=npoint;
else
M=N-(m-1)*tao;
end
Y=zeros(M,m);
for i=1:m
Y(:,i)=x((1:M)+(i-1)*tao)';
end
%

请先登录，再回答此问题。

回答（1 个）

SOUMNATH PAUL
SOUMNATH PAUL 2024-1-17
Hi,
Kindly do the following changes to your code to fix the issue:
  1. Remove ‘f’ in the lyarosenstein function because we are not going to use it further in the code.
  2. You can plot the raw divergence ‘d’ instead of the logarithm of the divergence ‘log(d)’.
  3. Include scripts to output the current verbose values which will help you in debugging.
Kindly go through the functions section of the below link for a related example:
https://www.mathworks.com/matlabcentral/fileexchange/38424-largest-lyapunov-exponent-with-rosenstein-s-algorithm
Hope it helps!
Regards,
Soumnath

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Computations 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by