MATLAB仿真：随参数b变化的最大Lyapunov指数图

MATLAB仿真：随参数b变化的最大Lyapunov指数图/n/n本文使用MATLAB仿真，基于四阶龙格库塔方法求解随机微分方程，利用Rosenstein方法计算最大Lyapunov指数，并绘制其随参数b变化的曲线图。该系统描述为：/n/n$$/dot{x} = /begin{bmatrix} x_2 // -sin(x_1) + /epsilon*(-ax_2 + c + bcos(kt)) + hcos(/omega_2 + /bar{/gamma}B(t)) /end{bmatrix}$$/n/n其中$$/xi(t)=h/cos(/omega_2+/bar{/gamma}B(t))$$, $$B(t)$$是单位维纳过程系统。/n/n参数设置：/n/n a = 0.02/n* c = 0.2/n* k = 1/n* h = 0.6/n* w2 = 1.8/n* gamma = 0.4/n* epsilon = 0.1/n/nMATLAB代码：/n/nmatlab/n% 参数设置/na = 0.02;/nc = 0.2;/nk = 1;/nh = 0.6;/nw2 = 1.8;/ngamma = 0.4;/nepsilon = 0.1;/n/n% 定义微分方程和噪声项/nf = @(t, x) [x(2); -sin(x(1)) + epsilon*(-a*x(2) + c + b*cos(k*t)) + h*cos(w2 + gamma*randn)];/ng = @(t, x) [0; 1];/n/n% 定义初始条件和计算时间/nx0 = [0.5; 0];/ntspan = [0 100];/nt = linspace(tspan(1), tspan(2), 10000);/n/n% 计算最大Lyapunov指数/nb_range = 0:0.05:2;/nlyapunov_max = zeros(size(b_range));/nfor i = 1:length(b_range)/n % 设置参数b/n b = b_range(i);/n /n % 计算轨道/n [t, x] = rk4_sde(f, g, x0, tspan, 0.01);/n /n % 计算雅可比矩阵和其特征值/n J = zeros(2, 2, length(t));/n lambda = zeros(2, length(t));/n for j = 1:length(t)/n J(:, :, j) = [0 1; -cos(x(j, 1)) epsilon*b*sin(k*t(j));];/n lambda(:, j) = eig(J(:, :, j));/n end/n /n % 计算最大Lyapunov指数/n lyapunov_max(i) = max(sum(log(abs(lambda)), 2)/(t(end) - t(1)));/nend/n/n% 画图/nfigure;/nplot(b_range, lyapunov_max);/nxlabel('Parameter b');/nylabel('Maximum Lyapunov exponent');/n/n/n其中，rk4_sde 是一个实现了四阶龙格库塔方法求解随机微分方程的函数，可以在以下链接中找到：https://www.mathworks.com/matlabcentral/fileexchange/71610-rk4_sde-m