MATLAB Program for Calculating Lyapunov Exponents

Below is the MATLAB source code for calculating Lyapunov exponents:

function [lyap, y_out] = lyapunov_exponent(x, f, df, n, d, a)

% x - Initial value (starting point for the system)

% f - System equation (function handle representing the dynamic system)

% df - Derivative of the system equation (Jacobian matrix or derivative function)

% n - Number of iterations (determines the accuracy and computation length)

% d - Step size (time step for numerical integration)

% a - System parameters (constants that define the specific system behavior)

y_out = zeros(1, n);

y_out(1) = log(norm(x));

for i = 2:n

y_out(i) = y_out(i-1) + d * df(x, a);

x = f(x, a);

end

lyap = y_out(n)/n;

end

This program implements the standard algorithm for computing Lyapunov exponents by iteratively tracking the exponential divergence rate of nearby trajectories in phase space. The core algorithm accumulates the logarithmic growth rate of perturbations over time, then averages the result to obtain the characteristic exponent. To use this program, you need to provide the initial condition, system dynamics function, its derivative (Jacobian), number of iterations, step size, and system parameters. The function returns both the final Lyapunov exponent and the complete evolution history, allowing you to analyze the system's chaotic properties and convergence behavior. A positive Lyapunov exponent typically indicates chaotic behavior in the system.