用MATLAB绘制混沌系统最大Lyapunov指数随参数b变化图

用MATLAB绘制混沌系统最大Lyapunov指数随参数b变化图

本文使用MATLAB编写代码，通过ODE算法和四阶龙格库塔法求解随机微分方程，绘制混沌系统最大Lyapunov指数随参数b变化的图像。该系统描述为一个非线性微分方程，包含随机噪声项，其最大Lyapunov指数能反映系统的混沌程度。

系统模型

该混沌系统由以下两个方程描述：

$\xi(t) = h*cos(\omega_2 + \bar{\gamma}*B(t))$

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

其中，

• $B(t)$ 是单位维纳过程

• $a=0.02, c=0.2, k=1, h=0.6, w_2=1.8, \bar{\gamma}=0.4, \epsilon=0.1$

MATLAB代码

% 定义随机微分方程的函数
function dx = chaotic_system(t,x,b)
a = 0.02;
c = 0.2;
k = 1;
h = 0.6;
w2 = 1.8;
gamma = 0.4;
epsilon = 0.1;

% 定义单位维纳过程
dB = randn();

% 计算xi(t)
xi = h*cos(w2 + gamma*dB);

% 计算dx/dt
dx = [x(2); -sin(x(1)) - epsilon*(-a*x(2) + c + b*cos(k*t)) + xi];
end

% 计算最大Lyapunov指数的函数
function lyap = lyapunov_exponent(b)
% 初始化
tspan = [0, 100];
x_0 = [0, 0.1*rand()];

% 设置ODE算法参数
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);

% 使用ODE算法求解随机微分方程
[t, x] = ode45(@(t,x) chaotic_system(t,x,b), tspan, x_0, options);

% 计算最大Lyapunov指数
n = length(t);
v = zeros(n, 2);
v(1, :) = [1, 0];
for i = 2:n
dx = chaotic_system(t(i), x(i, :), b);
J = [0, 1; -cos(x(i, 1)), -epsilon*(-a) - epsilon*b*k*sin(k*t(i))];
v(i, :) = J*v(i-1, :)'*norm(dx)/norm(v(i-1, :));
end
lyap = (1/n)*sum(log10(sqrt(sum(v.^2, 2))));
end

% 绘制随参数b变化的最大Lyapunov指数图
b_range = 0:0.1:10;
lyapunov = arrayfun(@lyapunov_exponent, b_range);
plot(b_range, lyapunov);
xlabel('b');
ylabel('Lyapunov exponent');
title('Lyapunov exponent vs b');

结果分析

运行该程序，可以得到如下Lyapunov指数图：

根据图像可以看出，当b的值在1.5左右时，系统的混沌程度最强，最大Lyapunov指数最大。当b的值继续增大时，系统的混沌程度逐渐减弱，最大Lyapunov指数也逐渐减小。