如何用MATLab计算有界噪声下的duffing振子微分方程含随机项并根据bennetin提出的算法画出其lyapunov指数谱图

以下是MATLab代码：

% 参数设定
a = 1;
b = 1;
gamma = 0.3;
omega = 1;
dt = 0.01;
t = 0:dt:100;
N = length(t);
x0 = [0.1 0];

% 随机项生成
sigma = 0.1;
dW = sigma*randn(2,N);

% 求解微分方程
x = zeros(2,N);
x(:,1) = x0;
for i = 2:N
    x1 = x(1,i-1);
    x2 = x(2,i-1);
    dW1 = dW(1,i-1);
    dW2 = dW(2,i-1);
    f = [x2; -a*x1-b*x1^3+gamma*cos(omega*t(i))];
    g = [0; 1];
    x(:,i) = x(:,i-1) + f*dt + g*dW(:,i-1);
end

% 计算lyapunov指数谱
m = 10;
L = zeros(1,m);
for j = 1:m
    delta = 1e-6;
    v0 = randn(2,1);
    v0 = v0/norm(v0);
    for i = 1:N
        x1 = x(1,i);
        x2 = x(2,i);
        f = [x2; -a*x1-b*x1^3+gamma*cos(omega*t(i))];
        g = [0; 1];
        df = [0 1; -3*b*x1^2-a 0];
        dg = [0; 0];
        v0 = v0 + (df*v0+dg*dW(:,i))*dt;
        v0 = v0/norm(v0);
        L(j) = L(j) + log(norm(v0+delta))/N;
    end
end

% 画出lyapunov指数谱图
plot(1:m, L, 'o-');
xlabel('指数');
ylabel('lyapunov指数');

代码中，首先设定了参数，包括振子的参数（a、b、gamma、omega）、时间步长（dt）、时间范围（t）、初始条件（x0）等。然后生成了随机项，采用标准正态分布的随机数乘以sigma。接着用欧拉方法求解微分方程，其中f是微分方程的右侧，g是随机项的系数。最后计算了lyapunov指数谱，采用了bennetin提出的算法，即在不同的初始方向上求解微分方程，计算每个方向上的lyapunov指数，最后求平均。最后画出了lyapunov指数谱图