其中$$xit=hcosomega_2+bargammaBt$$$$Bt$$是单位维纳过程系统$$ddotx+sinx-epsilon-adotx+c+bcoskt=xit$$取$$a=002c=02k=1h=06w_2=18bargamma=04epsilon=01$$用MATLAB根据Rosenstein 方法编写该随机系统的随参数b变化的最大lyapunov指数图最好是用四阶龙格库塔并将系统

以下是MATLAB代码：

function [t,lyapunov]=lyapunov_exponent()
% Parameters
a=0.02; c=0.2; k=1; h=0.6; w2=1.8; gamma=0.4; epsilon=0.1;
b_range=0:0.01:2; % range of b values
dt=0.01; % time step
t0=0; tf=1000; % initial and final time
tol=1e-6; % tolerance for ODE solver
% Initial conditions
x0=0; xp0=0; y0=0;
lyapunov=zeros(length(b_range),1);
for i=1:length(b_range)
    b=b_range(i);
    % Define the ODE system
    f=@(t,xp,x,y) [y+epsilon*a*xp-epsilon*c*x-epsilon*b*cos(k*t);-x+gamma*h*cos(w2+gamma*y);x];
    % Solve the ODE system using Runge-Kutta 4
    x=zeros(3,1);
    xp=zeros(3,1);
    y=zeros(3,1);
    x(1)=x0;
    xp(1)=xp0;
    y(1)=y0;
    for j=1:2
        t=t0+(j-1)*dt;
        k1=f(t,xp(j),x(j),y(j));
        k2=f(t+dt/2,xp(j)+dt*k1(1)/2,x(j)+dt*k1(2)/2,y(j)+dt*k1(3)/2);
        k3=f(t+dt/2,xp(j)+dt*k2(1)/2,x(j)+dt*k2(2)/2,y(j)+dt*k2(3)/2);
        k4=f(t+dt,xp(j)+dt*k3(1),x(j)+dt*k3(2),y(j)+dt*k3(3));
        xp(j+1)=xp(j)+(dt/6)*(k1(1)+2*k2(1)+2*k3(1)+k4(1));
        x(j+1)=x(j)+(dt/6)*(k1(2)+2*k2(2)+2*k3(2)+k4(2));
        y(j+1)=y(j)+(dt/6)*(k1(3)+2*k2(3)+2*k3(3)+k4(3));
    end
    % Calculate the tangent vector
    tangent_vector=[x(2)-x(1); xp(2)-xp(1); y(2)-y(1)];
    tangent_vector=tangent_vector/norm(tangent_vector);
    % Perturb the initial conditions for the next iteration
    x0=x(2)+tol*tangent_vector(2);
    xp0=xp(2)+tol*tangent_vector(1);
    y0=y(2)+tol*tangent_vector(3);
    % Calculate the Lyapunov exponent
    lyapunov(i)=1/(2*dt*i)*log(norm(tangent_vector)/tol);
end
% Plot the results
plot(b_range,lyapunov);
xlabel('b');
ylabel('Lyapunov exponent');
end

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