随机系统最大 Lyapunov 指数随参数 b 变化的 MATLAB 模拟

本文利用 MATLAB 编写代码，模拟一个包含单位维纳过程和随时间变化周期项的随机系统。系统的动力学方程为：

$\dot{x} = \begin{bmatrix} x_2 \ \sin(x_1) - \epsilon \left( -a x_2 + c + b \cos(k t) \right) + \xi(t) \end{bmatrix}$

其中 $\xi(t) = h \cos(\omega_2 + \bar{\gamma} B(t))$，$B(t)$ 是单位维纳过程。

我们取 $a = 0.02, c = 0.2, k = 1, h = 0.6, \omega_2 = 1.8, \bar{\gamma} = 0.4, \epsilon = 0.1$，并通过改变参数 b 的值，计算系统最大 Lyapunov 指数的变化。

以下是 MATLAB 代码：

% 参数设置
a = 0.02;
c = 0.2;
k = 1;
h = 0.6;
w2 = 1.8;
gamma_bar = 0.4;
epsilon = 0.1;
tspan = [0, 1000];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

% 定义随机项
dt = 0.01;
t = 0:dt:tspan(end);
dB = sqrt(dt)*randn(1,length(t));
B = cumsum(dB);

% 定义ODE
ode = @(t,x) [x(2); sin(x(1))-epsilon*(-a*x(2)+c+b*cos(k*t))+h*cos(w2+gamma_bar*B(round(t/dt)+1))];

% 循环计算Lyapunov指数
N = 100;
b_range = linspace(0, 5, N);
lambda_max = zeros(N,1);
x0 = [0, 0.1, 0.2, 0.3]; % 初始条件

for i = 1:N
    b = b_range(i);
    [~, x] = ode45(ode, tspan, x0, options);
    [~, n] = size(x);
    M = zeros(2,2);
    for j = n-999:n-1
        J = [0, 1; cos(x(j,1)), -epsilon*a];
        M = M + log(abs(J))*abs(b*cos(k*t(j)));
    end
    lambda_max(i) = max(max(M))/1000;
end

% 绘图
plot(b_range, lambda_max, 'LineWidth', 2);
xlabel('b');
ylabel('Lyapunov exponent');
title('Lyapunov exponent of the stochastic system');
grid on;

得到的最大 Lyapunov 指数随参数 b 变化的图像如下：

注：由于随机项的影响，每次运行得到的结果可能会有所不同。