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
You may refer to the MathWorks documentation for more information on “evalin” function. Here is the link to it:
I hope this helps!
