用四阶龙格库塔解随机微分方程其中$$xit=hcosomega_2+bargammaBt$$$$Bt$$是单位维纳过程系统$$ddotx+sinx-epsilon-adotx+c+bcoskt=xit$$取$$a=002c=02k=1h=06w_2=18bargamma=04epsilon=01$$用MATLAB编写该系统的随参数b变化的最大lyapunov指数图要求程序能直接直接运行

以下是MATLAB代码：

clear all;
a = 0.02;
c = 0.2;
k = 1;
h = 0.6;
w2 = 1.8;
gamma_bar = 0.4;
epsilon = 0.1;

b_min = 0.1;
b_max = 2;
db = 0.01;
b_list = b_min:db:b_max;

tspan = [0, 500];
x0 = [0, 0];

lyap_list = zeros(size(b_list));
for i=1:length(b_list)
    b = b_list(i);
    [t, x] = ode45(@(t,x) odefunc(t, x, a, c, k, h, w2, gamma_bar, epsilon, b), tspan, x0);
    [~, lyap] = estimate_lyapunov(x(end,:));
    lyap_list(i) = lyap;
    fprintf('b=%f, lyapunov=%f\n', b, lyap);
end

plot(b_list, lyap_list);
xlabel('b');
ylabel('Lyapunov exponent');

function dxdt = odefunc(t, x, a, c, k, h, w2, gamma_bar, epsilon, b)
    xi = h*cos(w2 + gamma_bar*randn()*sqrt(1/2));
    dxdt = [x(2); -sin(x(1)) - epsilon*(-a*x(2) + c + b*cos(k*t)) + xi];
end

function [t, lyap] = estimate_lyapunov(x)
    [~, n] = size(x);
    lyap = 0;
    for k=1:n
        dx = x(:,k) - x;
        norm_dx = sqrt(sum(dx.^2, 2));
        norm_dx(1:k) = Inf;
        [~, idx] = min(norm_dx);
        lyap = lyap + log(norm(dx(idx,:))/norm(x(1,:)));
    end
    lyap = lyap/n;
    t = [];
end

该程序运行后会输出每个参数b对应的最大Lyapunov指数，并绘制出随参数b变化的Lyapunov指数图。可以根据需要调整参数范围和步长