how to find lyapunov exponent of a 3d discrete system in matlab?

for example like predator-prey model

Answer 1

[Divyam]

Hi @ik kumar,

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.