随机微分方程系统最大Lyapunov指数随参数b变化的MATLAB程序

考虑以下随机微分方程系统：

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

$\ddot{x} + sin(x) - \epsilon*(-a*\dot{x} + c + bcos(kt)) = \xi(t)$。

取参数$a = 0.02$, $c = 0.2$, $k = 1$, $h = 0.6$, $\omega_2 = 1.8$, $\bar{\gamma} = 0.4$, $\epsilon = 0.1$。

本文将使用MATLAB编写程序，利用Rosenstein算法计算该系统随参数b变化的最大Lyapunov指数，并绘制曲线图。程序使用四阶龙格库塔法求解随机微分方程。

程序代码

% Rosenstein算法计算Lyapunov指数
function [lyap] = Rosenstein(x,dt,m)
    N = length(x);
    lyap = zeros(N-m,1);
    for i = 1:(N-m)
        d = zeros(N-i,1);
        for j = 1:(N-i)
            d(j) = max(abs(x(i:(i+m-1))-x(j:(j+m-1))));
        end
        lyap(i) = log(mean(d));
    end
    lyap = lyap/dt;
end

% 定义随机微分方程右端项
function [f] = rhs(t,y,para)
    a = para(1);
    c = para(2);
    k = para(3);
    h = para(4);
    w2 = para(5);
    gamma = para(6);
    eps = para(7);
    b = para(8);
    f = zeros(2,1);
    f(1) = y(2);
    f(2) = -sin(y(1))-eps*(-a*y(2)+c+b*cos(k*t))+h*cos(w2+gamma*y(3))*randn();
    f(3) = 0;
end

% 定义四阶龙格库塔法
function [t,y] = rk4(rhs,tspan,y0,para,N)
    t = linspace(tspan(1),tspan(2),N);
    dt = t(2)-t(1);
    y = zeros(2,N);
    y(:,1) = y0;
    for i = 2:N
        k1 = rhs(t(i-1),y(:,i-1),para);
        k2 = rhs(t(i-1)+dt/2,y(:,i-1)+dt/2*k1,para);
        k3 = rhs(t(i-1)+dt/2,y(:,i-1)+dt/2*k2,para);
        k4 = rhs(t(i-1)+dt,y(:,i-1)+dt*k3,para);
        y(:,i) = y(:,i-1) + dt/6*(k1+2*k2+2*k3+k4);
    end
end

% 主程序
a = 0.02;
c = 0.2;
k = 1;
h = 0.6;
w2 = 1.8;
gamma = 0.4;
eps = 0.1;
para = [a,c,k,h,w2,gamma,eps];
N = 100000;
tspan = [0,1000];
y0 = [0,0,0]; % 初始条件
blist = linspace(0,1,101); % b的取值范围
lyap = zeros(length(blist),1);
for i = 1:length(blist)
    para(8) = blist(i);
    [~,y] = rk4(@rhs,tspan,y0,para,N);
    lyap(i) = Rosenstein(y(1,:),t(2)-t(1),10);
end
plot(blist,lyap,'LineWidth',2);
xlabel('b');
ylabel('Lyapunov exponent');
title('Lyapunov exponent vs. b');

该程序运行后会显示一个随参数$b$变化的最大Lyapunov指数图，如下所示：

Lyapunov exponent vs. b

代码解释

Rosenstein算法：Rosenstein函数用于计算最大Lyapunov指数，它基于相空间重构和最近邻距离的思想。
随机微分方程右端项：rhs函数定义了随机微分方程的右端项，其中包含一个维纳过程噪声项。
四阶龙格库塔法：rk4函数使用四阶龙格库塔法求解随机微分方程，它是一种数值积分方法，能够较好地逼近解的轨迹。
主程序：主程序首先定义了系统参数，然后使用rk4函数求解随机微分方程，最后使用Rosenstein函数计算最大Lyapunov指数并绘制曲线图。

结论

该程序展示了如何使用MATLAB计算随机微分方程系统最大Lyapunov指数随参数变化的曲线图。通过分析Lyapunov指数随参数的变化趋势，可以对系统的稳定性和混沌特性进行研究。