How I cant Plot the Spectrum of Lyapunov exponents of system versus the parameter c∈ [0, 200], with a = 10, b = 28 ?

Question

YOUSSEF El MOUSSATI 2024-4-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2109156-how-i-cant-plot-the-spectrum-of-lyapunov-exponents-of-system-versus-the-parameter-c-0-200-with

评论： YOUSSEF El MOUSSATI 2024-9-5，12:30

Greetings, all. I trust you're doing well. My intention is to generate a spectrum depicting the Lyapunov exponents of the system concerning the parameter c within the range of \( [0, 200] \), given ( a = 10 ) and ( b = 28 ). Fig(11.png).The system equations are as follows:

x'= -a * x + b * y - y * z;

y' = x + z * x;

z' = -c * z + y^2;

I'm employing specific codes to calculate the Lyapunov exponents for c = 2 Fig(22.png).

codes :

%%%%%%%%% system_equations %%%%%%%%%%%%%%%%%%%
% Define the system equations
function f = system_equations(t, X)
% Define the parameters
a = 10;
b = 28;
c = 2;
% Extract variables
x = X(1);
y = X(2);
z = X(3);
% Extract the Jacobian matrix Y from the extended state vector X
Y = [X(4), X(7), X(10);
X(5), X(8), X(11);
X(6), X(9), X(12)];
% Initialize the output vector
f = zeros(9, 1);
% System equations
f(1) = -a * x + b * y - y * z;
f(2) = x + z * x;
f(3) = -c * z + y^2;
% Jacobian matrix
Jac = [-a, b - z, -y;
1 + z, 0, x;
0, 2 * y, -c];
% Variational equation
f(4:12) = Jac*Y;
end

%%%%%%%%% systemequations_Lyapunov %%%%%%%%%%%%%%%%%%%
function [Texp, Lexp] = systemequations_Lyapunov(n, rhs_ext_fcn, fcn_integrator, tstart, stept, tend, ystart, ioutp)
clc;
close all;
n = 3;
rhs_ext_fcn = @system_equations;
fcn_integrator = @ode45;
tstart = 0;
stept = 0.1;
tend = 1000;
ystart = [7.2 7.8 2.3];
ioutp = 10;
n1 = n;
n2 = n1*(n1 + 1);
% Number of steps.
nits = round((tend - tstart)/stept);
% Memory allocation.
y = zeros(n2, 1);
cum = zeros(n1, 1);
y0 = y;
gsc = cum;
znorm = cum;
% Initial values.
y(1:n) = ystart(:);
for i = 1:n1
y((n1 + 1)*i) = 1.0;
end
t = tstart;
% Main loop.
for iterlya = 1:nits
% Solutuion of extended ODE system.
[T, Y] = feval(fcn_integrator, rhs_ext_fcn, [t t + stept], y);
t = t + stept;
y = Y(size(Y, 1),:);
for i = 1:n1
for j = 1:n1
y0(n1*i + j) = y(n1*j + i);
end
end
% Construct new orthonormal basis by Gram-Schmidt.
znorm(1) = 0.0;
for j = 1:n1
znorm(1) = znorm(1) + y0(n1*j+1)^2;
end
znorm(1) = sqrt(znorm(1));
for j = 1:n1
y0(n1*j + 1) = y0(n1*j + 1)/znorm(1);
end
for j = 2:n1
for k = 1:(j - 1)
gsc(k) = 0.0;
for l = 1:n1
gsc(k) = gsc(k) + y0(n1*l + j)*y0(n1*l + k);
end
end
for k = 1:n1
for l = 1:(j - 1)
y0(n1*k + j) = y0(n1*k + j) - gsc(l)*y0(n1*k + l);
end
end
znorm(j) = 0.0;
for k = 1:n1
znorm(j) = znorm(j) + y0(n1*k+j)^2;
end
znorm(j)=sqrt(znorm(j));
for k = 1:n1
y0(n1*k + j) = y0(n1*k + j)/znorm(j);
end
end
% Update running vector magnitudes.
for k = 1:n1
cum(k) = cum(k) + log(znorm(k));
end
% Normalize exponent.
for k = 1:n1
lp(k) = cum(k)/(t - tstart);
end
% Output modification.
if iterlya == 1
Lexp = lp;
Texp = t;
else
Lexp = [Lexp; lp];
Texp = [Texp; t];
end
for i = 1:n1
for j = 1:n1
y(n1*j + i) = y0(n1*i + j);
end
end
end
% Show the Lyapunov exponent values on the graph.
str1 = num2str(Lexp(nits, 1));
str2 = num2str(Lexp(nits, 2));
str3 = num2str(Lexp(nits, 3));
pl = plot(Texp, Lexp,'LineWidth',2);
pl(1).Color = 'b';
pl(2).Color = 'g';
pl(3).Color = 'r';
legend('\lambda_1', '\lambda_2','\lambda_3')
legend('Location','northeast')
title('Dynamics of Lyapunov Exponents');
text(205, 1.5, '\lambda_1 = ','Fontsize',20);
text(220, 1.68, str1,'Fontsize',20);
text(205, -1, '\lambda_2 = ','Fontsize',20);
text(220, -0.73, str2,'Fontsize',20);
text(205, -13.8, '\lambda_3 = ','Fontsize',20);
text(220, -13.5, str3,'Fontsize',20);
xlabel('Time');
ylabel('Lyapunov Exponents');
axis([-1,300, -16,6]);
set(gca,'FontSize',20)
end
% End of plot

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Aastha 2024-9-3，2:54

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2109156-how-i-cant-plot-the-spectrum-of-lyapunov-exponents-of-system-versus-the-parameter-c-0-200-with#answer_1509339

在 MATLAB Online 中打开

Hi YOUSSEF El MOUSSATI,

I understand that you want to generate a plot of Lyapunov exponents as a function of the parameter “c” given a = 10 and b = 28. You can follow the steps mentioned below to plot Lyapunov exponents as a function of parameter “c”:

Define the range for parameter “c”: Define the range over which the parameter “c” will vary. Based on the graph attached, the range of “c” is from 0 to 200 with a step size of 20.
Compute Lyapunov Exponents: You can use a for loop to iterate over the range of values of “c”. For each value of “c”, compute and save the three Lyapunov exponents for each time step.
Modify the Solver: To vary “c” inside the “ode45” solver without changing the “system_equations” function definition, you can use the MATLAB function “evalin”.
Inside the “system_equations” function, you can set the value of parameter “c” from the base workspace to the current value of “c” that is being used in the for loop.
Plot Lyapunov Exponents versus “c”: After computing the Lyapunov exponents as a function of time steps for all “c” values, you can select a specific time step and use the corresponding Lyapunov exponent values to generate the plot of Lyapunov exponents versus “c”
You can experiment with different time steps to match the attached figure as closely as possible.

I am hereby attaching a sample code to plot Lyapunov exponent versus “c” for a = 10 and b = 28 following the above-mentioned steps:

close all
clear
clc
crange=0:20:200;
Lyapunov_exp = zeros(length(crange),10000,3);
cnt = 1;
for c = crange
n = 3;
rhs_ext_fcn = @system_equations;
fcn_integrator = @ode45;
tstart = 0;
stept = 0.1;
tend = 1000;
ystart = [7.2 7.8 2.3];
ioutp = 10;
n1 = n;
n2 = n1*(n1 + 1);
% Number of steps.
nits = round((tend - tstart)/stept);
% Memory allocation.
y = zeros(n2, 1);
cum = zeros(n1, 1);
y0 = y;
gsc = cum;
znorm = cum;
% Initial values.
y(1:n) = ystart(:);
for i = 1:n1
y((n1 + 1)*i) = 1.0;
end
t = tstart;
% Main loop.
for iterlya = 1:nits
% Solutuion of extended ODE system.
[T, Y] = feval(fcn_integrator, rhs_ext_fcn, [t t + stept], y);
t = t + stept;
y = Y(size(Y, 1),:);
for i = 1:n1
for j = 1:n1
y0(n1*i + j) = y(n1*j + i);
end
end
% Construct new orthonormal basis by Gram-Schmidt.
znorm(1) = 0.0;
for j = 1:n1
znorm(1) = znorm(1) + y0(n1*j+1)^2;
end
znorm(1) = sqrt(znorm(1));
for j = 1:n1
y0(n1*j + 1) = y0(n1*j + 1)/znorm(1);
end
for j = 2:n1
for k = 1:(j - 1)
gsc(k) = 0.0;
for l = 1:n1
gsc(k) = gsc(k) + y0(n1*l + j)*y0(n1*l + k);
end
end
for k = 1:n1
for l = 1:(j - 1)
y0(n1*k + j) = y0(n1*k + j) - gsc(l)*y0(n1*k + l);
end
end
znorm(j) = 0.0;
for k = 1:n1
znorm(j) = znorm(j) + y0(n1*k+j)^2;
end
znorm(j)=sqrt(znorm(j));
for k = 1:n1
y0(n1*k + j) = y0(n1*k + j)/znorm(j);
end
end
% Update running vector magnitudes.
for k = 1:n1
cum(k) = cum(k) + log(znorm(k));
end
% Normalize exponent.
for k = 1:n1
lp(k) = cum(k)/(t - tstart);
end
% Output modification.
if iterlya == 1
Lexp = lp;
Texp = t;
else
Lexp = [Lexp; lp];
Texp = [Texp; t];
end
for i = 1:n1
for j = 1:n1
y(n1*j + i) = y0(n1*i + j);
end
end
end
Lyapunov_exp(cnt,:,:)=Lexp;
cnt = cnt + 1;
end
% Plot the variation as a function of c
timeStep = 10; % choose a time step to plot the variation
plotLexp = zeros(length(crange),3);
for cnt = 1:length(crange)
    Lexp = Lyapunov_exp(cnt, timeStep, :);
    plotLexp(cnt,:)=squeeze(Lexp);
end
figure(1);
for lambda_cnt = 1:3
    plot(crange, plotLexp(:,lambda_cnt),'DisplayName',['lambda ',num2str(cnt)]);
    hold on
end
hold off
grid on
xlabel('c');
ylabel('Lyapunov Exponents');
legend
title('Variation of Lyapunov Exponents v.s c');

% Define the system equations
function f = system_equations(t, X)
% Define the parameters
a = 10;
b = 28;
c=evalin('base','c');
% Extract variables
x = X(1);
y = X(2);
z = X(3);
% Extract the Jacobian matrix Y from the extended state vector X
Y = [X(4), X(7), X(10);
X(5), X(8), X(11);
X(6), X(9), X(12)];
% Initialize the output vector
f = zeros(9, 1);
% System equations
f(1) = -a * x + b * y - y * z;
f(2) = x + z * x;
f(3) = -c * z + y^2;
% Jacobian matrix
Jac = [-a, b - z, -y;
1 + z, 0, x;
0, 2 * y, -c];
% Variational equation
f(4:12) = Jac*Y;
end