Dear Mr Elhadi

Here is an example of MATLAB code to calculate the Lyapunov exponent (LE) of a discrete chaotic system in 2D:

```matlab

% Define the system equations

function dxdt = chaoticSystem(t, x)

a = 1.4; % Parameter value

b = 0.3; % Parameter value

dxdt = [a * (x(2) - x(1));

x(1) * (b - x(2)) - x(1) * x(3);

x(1) * x(2) - b * x(3)];

end

% Numerically solve the system of equations using ode45

tspan = [0 100]; % Time span for integration

x0 = [1; 1; 1]; % Initial conditions

options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6); % Solver options

[t, x] = ode45(@chaoticSystem, tspan, x0, options);

% Estimate the Lyapunov exponent

m = length(t); % Number of time steps

n = length(x0); % Dimension of the system

LEs = zeros(1, n); % Initialize Lyapunov exponents

delta0 = 1e-6; % Initial perturbation magnitude

perturbation = eye(n) * delta0;

for i = 1:m

J = jacobianMatrix(x(i, :)); % Calculate Jacobian matrix

perturbation = J * perturbation; % Propagate perturbation

[~, R] = qr(perturbation); % QR factorization to orthogonalize perturbation

perturbation = R; % Update perturbation

norm_perturbation = norm(perturbation, 'fro'); % Frobenius norm of perturbation

perturbation = perturbation / norm_perturbation; % Normalize perturbation

LEs = LEs + log(norm_perturbation) / (t(i) - t(1)); % Update Lyapunov exponents

end

% Display Lyapunov exponents

disp("Lyapunov Exponents:");

disp(LEs);

% Function to compute the Jacobian matrix of a given state

function J = jacobianMatrix(x)

a = 1.4; % Parameter value

b = 0.3; % Parameter value

J = [-a, a, 0;

b - x(2) - x(3), -1, -x(1);

x(2), x(1), -b];

end

```

In this code, the chaotic system is defined by the `chaoticSystem` function, which represents the system of differential equations. The `ode45` solver is used to numerically integrate the system over a given time span.

After obtaining the solution trajectory, the Lyapunov exponents are estimated using the tangent space method. The Jacobian matrix of the system is computed at each time step using the `jacobianMatrix` function. Perturbations are propagated along the trajectory, orthogonalized using QR factorization, and normalized to maintain stability. The Lyapunov exponents are updated based on the rate of perturbation growth.

Finally, the Lyapunov exponents are displayed as output using the `disp` function. The number of Lyapunov exponents corresponds to the dimensionality of the system.

And in the link below, you may find your answer

https://www.google.com/url?sa=t&source=web&rct=j&opi=89978449&url=https://www.researchgate.net/post/How-to-compute-both-LEs-for-2D-discrete-chaotic-logistic-map-in-MATLAB&ved=2ahUKEwjOuMHa_ZqCAxWTtqQKHZdcAC0QFnoECCcQAQ&usg=AOvVaw14-0LfJtvtxUgd5vyG4n6I

One thing to add, hope it helps in this case its might,

I believe the default type is floating point in the matlab code and the ODE45 is slow.

Article Lyapunov Exponent

as mentioned above in wiki

for hard constraints, stability and determinism, reversibility using fixed point, if float is the default. for reproducibility if it can be modified simply.

here is the official way.

https://www.mathworks.com/help/fixedpoint/ug/create-fixed-point-data.html

this link below you can see that even a very crude integration will be stable and reversible , you can see how quickly a 2 hinge pendulum goes into chaos caused just by the floating point errors , arithmetic non associativity and such.

and how he convert binary int representation to fixed over a domain: without the fixed types, by dividing an unsigned int length by a read number for the domain.

https://www.josstam.com/reversible see even the rough integrators are stable and deterministic

Hi Mehallel Elhadi

this matlab code i hope so help you:

function LEs = henonLEs(a, b, N)

% a and b are parameters of the Henon map

% N is the number of iterations

% Initial conditions

x = 0.1;

y = 0.1;

% Initial vectors for the tangent map

v1 = [1;0];

v2 = [0;1];

sum_lambda1 = 0;

sum_lambda2 = 0;

for n = 1:N

% Compute the Jacobian matrix at (x,y)

J = [ -2*a*x, 1; b, 0];

% Evolve the vectors using the tangent map

v1_new = J * v1;

v2_new = J * v2;

% Gram-Schmidt orthogonalization

v1_new = v1_new / norm(v1_new);

v2_new = v2_new - (v2_new' * v1_new) * v1_new;

v2_new = v2_new / norm(v2_new);

% Compute the growth rate for the vectors

sum_lambda1 = sum_lambda1 + log(norm(v1_new));

sum_lambda2 = sum_lambda2 + log(norm(v2_new));

% Update the vectors for the next iteration

v1 = v1_new;

v2 = v2_new;

% Evolve the system

x_new = y + 1 - a * x^2;

y = b * x;

x = x_new;

end

lambda1 = sum_lambda1 / N;

lambda2 = sum_lambda2 / N;

LEs = [lambda1, lambda2];

end

see this

https://www.mathworks.com/help/predmaint/ref/lyapunovexponent.html

This can be readily adjusted for 2D recurrences

https://github.com/sahariml/Lyapunov_diag_bif_3D-PLM-Map

Thank you so much