当前位置：首页 > news >正文

基于预估校算法的分数阶混沌系统Lyapunov指数计算

news 2026/3/27 4:02:55

基于预估校算法（Adams-Bashforth-Moulton方法）的分数阶混沌系统Lyapunov指数计算

一、算法原理框架

数学模型：

采用Grünwald-Letnikov定义：

二、Matlab核心代码实现

%% 分数阶Lorenz系统参数设置
alpha = 0.95; % 分数阶阶数
sigma = 10; rho = 28; beta = 8/3;
h = 0.01; T = 1000; N = T/h;%% 预估校算法初始化
X = zeros(3,N); X(:,1) = [1;1;1]*0.1;
Y = zeros(3,N); Y(:,1) = X(:,1) + 1e-6*randn(3,1);%% Lyapunov指数计算初始化
LE = zeros(3,1); LE(1) = 0.01; % 初始猜测
P = eye(3); % 线性化矩阵%% 主循环
for n = 1:N-1% 预估步 (Adams-Bashforth)f1 = fLorenz(X(:,n), sigma, rho, beta);f2 = fLorenz(X(:,n-1), sigma, rho, beta);X_pred = X(:,n) + h^alpha/2 * (3*f1 - f2);% 校正步 (Adams-Moulton)f_pred = fLorenz(X_pred, sigma, rho, beta);X(:,n+1) = X(:,n) + h^alpha/2 * (f1 + 3*f_pred);% 扰动演化Y(:,n+1) = X(:,n+1) + P * (Y(:,n) + 1e-6*randn(3,1));% 线性化更新P = P + h^alpha * (P*fLorenz(X(:,n+1), sigma, rho, beta)' ...- fLorenz(X(:,n+1), sigma, rho, beta)*P);% Lyapunov指数更新LE = LE + log(norm(Y(:,n+1))/norm(Y(:,n)));
end%% 结果归一化
LE = LE / T;%% 分数阶Lorenz系统定义
function dx = fLorenz(x, sigma, rho, beta)dx = zeros(3,1);dx(1) = sigma*(x(2)-x(1))^alpha;dx(2) = x(1)*(rho-x(3))^alpha - x(2);dx(3) = x(1)*x(2) - beta*x(3)^alpha;
end

三、关键算法改进

1. 分数阶导数计算优化

% 改进的Grünwald-Letnikov离散化
function d = gl_deriv(x, alpha)n = length(x);d = zeros(n,1);for k = 1:n-1d(k) = (d(k-1)*(1+alpha-k) + x(n-k+1)) / gamma(alpha+1);endd(n) = 0;
end

2. 扰动向量正交化

% Gram-Schmidt正交化过程
function Y = gram_schmidt(Y)[m,n] = size(Y);for i = 2:nY(:,i) = Y(:,i) - Y(:,1:i-1)*Y(:,1:i-1)'*Y(:,i);Y(:,i) = Y(:,i)/norm(Y(:,i));end
end

四、可视化分析

%% 三维相图绘制
figure;
plot3(X(1,:), X(2,:), X(3,:), 'b.');
hold on;
plot3(Y(1,:), Y(2,:), Y(3,:), 'r.');
xlabel('x'); ylabel('y'); zlabel('z');
title('分数阶Lorenz系统相空间轨迹');%% Lyapunov指数谱
figure;
bar(LE);
set(gca,'XTickLabel',{'Lyap1','Lyap2','Lyap3'});
xlabel('Lyapunov指数序号');
ylabel('指数值');
title('分数阶Lorenz系统Lyapunov指数谱');

五、应用扩展

1. 超混沌系统分析

% 4D分数阶Chen系统
function dx = fChen(x, a,b,c,d)dx = zeros(4,1);dx(1) = a*(x(2)-x(1))^alpha;dx(2) = (c-a)*x(1) - x(1)*x(3) + c*x(2) + d*x(4);dx(3) = x(2)*x(1) - b*x(3)^(2-alpha);dx(4) = (d+b)*x(4) - x(3)*x(2);
end

2. 动态分岔分析

% 分岔参数扫描
alpha_values = 0.85:0.01:0.99;
for i = 1:length(alpha_values)alpha = alpha_values(i);[~, LE] = compute_lyapunov(X0, alpha); % 调用主函数plot(alpha, LE(1),'bo'); hold on;
end
xlabel('分数阶阶数α'); ylabel('最大Lyapunov指数');
title('分数阶Lorenz系统分岔图');