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

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

Answer 1

[Aastha]

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”:

1. 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.
2. 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.
3. Modify the Solver: To vary “c” inside the “ode45” solver without changing the “system_equations” function definition, you can use the MATLAB function “evalin”.
4. 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.
5. 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”
6. 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

You may refer to the MathWorks documentation for more information on “evalin” function. Here is the link to it:

https://www.mathworks.com/help/simulink/slref/simulink.modelworkspace.evalin.html

I hope this helps!

[YOUSSEF El MOUSSATI]

Thank you 🙏