$$xit=hcosomega_2+bargammaBt$$$$Bt$$是单位维纳过程系统$$ddotx+sinx-epsilon-adotx+c+bcoskt=xit$$取$$a=002c=02k=1h=06w_2=18bargamma=04epsilon=01$$用MATLAB编写该系统的随参数b变化的最大lyapunov指数图用Rosenstein算法要求程序能直接直接运行

以下是MATLAB代码：

function lyapunov_exponent()
    a = 0.02;
    c = 0.2;
    k = 1;
    h = 0.6;
    w2 = 1.8;
    gamma_bar = 0.4;
    eps = 0.1;
    tspan = [0, 100];
    x0 = [0, 0];
    b_range = linspace(0.1, 2, 100); % b的范围
    lambda_max = zeros(size(b_range)); % 存储最大Lyapunov指数
    options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
    for i = 1:length(b_range)
        b = b_range(i);
        [~, ~, lambda] = ode45(@(t, x) odefun(t, x, a, c, k, h, w2, gamma_bar, eps, b), tspan, x0, options);
        lambda_max(i) = max(lambda);
    end
    plot(b_range, lambda_max);
    xlabel('b');
    ylabel('Lyapunov exponent');
end

function dxdt = odefun(t, x, a, c, k, h, w2, gamma_bar, eps, b)
    B = sqrt(2) * randn(); % 单位维纳过程
    xi = h * cos(w2 + gamma_bar * B);
    dxdt = [x(2); eps * (-a * x(2) + c + b * cos(k * t)) - sin(x(1)) + xi];
end

运行该函数即可得到随参数b变化的最大Lyapunov指数图