To compute the Lyapunov exponents of a 3D discrete dynamical system you would need to follow a numerical procedure on an evolving Jacobian matrix along a trajectory while performing QR decomposition at each step.
The following code illustrates the approach mentioned above:
function lyap3d
% Parameters (example values)
r = 1.1; a = 0.1; b = 0.2;
N = 10000; % iterations
dt = 1; % discrete system
% Initial condition
x = [0.5; 0.4; 0.3];
% Initialize orthonormal vectors
Q = eye(3);
le = zeros(1, 3); % To store Lyapunov exponents
for i = 1:N
% Evolve system
x = f(x, r, a, b);
% Compute Jacobian at current x
J = jacobian_f(x, r, a, b);
% Tangent map
Z = J * Q;
% QR decomposition
[Q, R] = qr(Z);
% Accumulate log of absolute values of diagonal
le = le + log(abs(diag(R)))';
end
% Average to get the Lyapunov exponents
le = le / (N * dt);
fprintf('Lyapunov Exponents:\n %f\n %f\n %f\n', le);
end
function xn = f(x, r, a, b)
% Example map function
xn = [
x(1)*(1 + r - a*x(2));
x(2)*(1 - b + a*x(1) - 0.1*x(3));
x(3)*(1 - 0.05*x(1) + 0.02*x(2));
];
end
function J = jacobian_f(x, r, a, b)
% Jacobian matrix of the map f
J = zeros(3);
J(1,1) = 1 + r - a*x(2);
J(1,2) = -a*x(1);
J(1,3) = 0;
J(2,1) = a*x(2);
J(2,2) = 1 - b + a*x(1) - 0.1*x(3);
J(2,3) = -0.1*x(2);
J(3,1) = -0.05*x(3);
J(3,2) = 0.02*x(3);
J(3,3) = 1 - 0.05*x(1) + 0.02*x(2);
end
Note: You can define the map function (f(x)) and the Jacobian matrix of partial derivatives (J(x)) as per your requirements.
