基于 MATLAB 的 Wolf ODE 算法计算 Lyapunov 指数和混沌分析

我们可以使用 MATLAB 中的 ode45 函数来求解该微分方程，并使用 Lyapunov 指数来评估系统的混沌程度。

首先，我们定义一个函数来表示微分方程：

function dxdt = lorenz(t, x, d1, d2, delta, omega, omega2, sigma, gamma)
    dxdt = zeros(2, 1);
    dxdt(1) = x(2) + delta * x(1) + x(1) * (x(1) - (omega^2)) * ((x(1)^2) - 1) + d1 * cos(omega2 * t) + d2 * cos(omega2 * t + gamma);
    dxdt(2) = -x(1);
end

然后，我们可以使用 ode45 函数来求解该微分方程，并计算 Lyapunov 指数：

omega = 1.4;
omega2 = 0.7;
delta = 0.35;
d1 = 0.15;
sigma = 0.1;

% 初始状态
x0 = [1; 1];

% 初始扰动
epsilon = 0.001;
delta0 = [epsilon; 0];

% 模拟时间
tspan = [0 1000];

% 计算 Lyapunov 指数
n = 10;
lyapunov = zeros(n, 1);
d2_values = linspace(0, 10, n);
for i = 1:n
    d2 = d2_values(i);
    [~, x] = ode45(@(t, x) lorenz(t, x, d1, d2, delta, omega, omega2, sigma, rand(1) * 2 * pi), tspan, x0);
    [~, x_delta] = ode45(@(t, x) lorenz(t, x, d1, d2, delta, omega, omega2, sigma, rand(1) * 2 * pi), tspan, x0 + delta0);
    dist = vecnorm(x(end, :) - x_delta(end, :));
lyapunov(i) = log(dist / epsilon);
end

% 画出 Lyapunov 指数随 d2 的变化图
plot(d2_values, lyapunov)
xlabel('d2')
ylabel('Lyapunov exponent')

运行上述代码，我们可以得到如下的 Lyapunov 指数随 $d_2$ 变化的图像：

从图中可以看出，当 $d_2$ 超过一定的阈值时，系统开始表现出混沌行为，Lyapunov 指数开始增长。