用于确定分数阶系统（FOS）的Lyapunov指数谱，包括分数阶Lorenz系统、4D分数阶Chen系统和分数阶Duffing振荡器（Matlab代码实现）

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

💥1 概述

摘要
Lyapunov指数提供了确定动力系统稳定性和分类极限集的定量证据。目前有几种成熟的技术可以计算整数阶系统的Lyapunov指数，然而，由于分数阶导数的非局部性，这些技术无法推广到分数阶系统。本文提出了一种用于确定分数阶系统Lyapunov指数谱的方法。所提出的方法是根据Grünwald-Letnikov导数的记忆原理严格推导而来的，因此它通常适用，甚至与整数阶系统兼容。分数阶Lorenz系统、分数阶Duffing振荡器和4维分数阶Chen系统等三个经典示例分别用于展示所提方法对不可测、非自治和低有效阶系统以及超混沌系统的有效性。模拟结果表明，所提方法在准确性和正确性上确实优于现有方法。

1. 引言
Lyapunov指数（LE）是由Oseledets [1]在他的多重遍历定理中引入的。基于Oseledets的理论，Benettin等人[2]首次报道了计算动力系统所有LE的方法。后来，Wolf等人[3]改进了Benettin的方法，并首次提出了基于Takens的重构技术[4]的时间序列估计LE的方法，后者在实验研究中被广泛使用。除了Benettin-Wolf算法外，过去四十年来已经报道了各种LE计算方案，主要包括从控制方程中确定[5]，[6]和从时间序列中估计[7]，[8]。此外，研究人员还开发了扰动方法，如扰动向量方法[9]和克隆动力学方法[10]，以及同步方法[11]。这些方法避免直接计算雅可比矩阵或求解变分方程，因此它们可以应用于更多情况，例如一般具有病态雅可比矩阵的非光滑系统。目前，上述方法已经有效应用于几个领域的动态特征化[12]，[13]，[14]，[15]。分数微积分（FC）是微积分理论的一个古老分支。在早期，与古典微积分不同，FC在纯数学理论上有更多的研究，但在物理学中只有零星的应用（由Valério等人总结[16]）。在过去的二十年里，更多的研究人员开始关注FC的跨学科应用，因为一些问题可以通过FC以新颖的方式进行优雅建模，例如粘弹性材料的机械行为[17]，[18]，[19]，图像加密[20]，[21]和流行病建模[22]，[23]。最近，Diethelm等人[24]和Sun等人[25]分别回顾了FC在物理学和工程学中的应用。

一般来说，分数阶微分方程的解不能像半群那样定义动力系统[26]，然而，这并不意味着分数阶微分方程与其相流之间的关系不能像整数阶微分方程一样建立起来。作为相流收敛或发散速率的度量，LE仍然是研究分数阶系统（FOS）动态演化的强大数学工具。Li等人[27]首次引入了FOS的LE的严格数学定义，并提出了确定其上下界的方法。此外，还开发了一些新方法，如时间序列方法[28]，[29]和扩展的Benettin-Wolf算法[30]，用于估计FOS的LE。在这些进展中，确认LE仍然有效用于确定稳定性、分形维数和极限集。然而，存在一些限制，例如，扩展的Benettin-Wolf算法更适用于准整数阶系统，因为非局部性不够显著而可以忽略。事实上，FOS的非局部性意味着相流的收敛和/或发散是历史相关的，因此在FOS中很难推广整数阶系统中的成熟方法。因此，通过将所有维度中收敛（或发散）的历史相关关系形式化，可以正确确定FOS的所有LE，这是本工作的主要贡献。详细文章见第4部分。

这篇文章中的MATLAB代码用于确定分数阶系统（FOS）的Lyapunov指数谱，包括三个经典示例：分数阶Lorenz系统、4D分数阶Chen系统和分数阶Duffing振荡器。该算法基于分数阶导数的记忆原理，对系统的维度和阶数没有限制。当阶数设置为1时，数值方法自动简化为前向Euler方案，因此该程序也可用于确定整数阶系统的Lyapunov指数。

亮点：
- 通用适用性，并与标准的整数阶系统兼容。
- 更好的准确性和正确性。
- 无需手动调整。
- 适用于不可测、非自治和低有效阶系统以及超混沌系统。

用于确定分数阶系统（FOS）的Lyapunov指数谱的研究综述

一、研究背景与意义

分数阶系统（FOS）通过引入分数阶微积分，能够更精确地描述具有记忆效应和长程相关性的复杂系统，如粘弹性材料、电力系统、生物神经网络等。Lyapunov指数谱是衡量动力系统稳定性和混沌特性的核心工具，其正负和数值大小直接反映系统对初值的敏感性、轨道分离率及吸引子结构。然而，分数阶导数的非局部性和奇异性使得传统整数阶系统的Lyapunov指数计算方法无法直接推广，需结合分数阶微积分的特性进行改进。

二、分数阶系统Lyapunov指数谱的计算方法

当前研究主要围绕以下两类方法展开，结合分数阶系统的特性进行优化：

1. 基于QR分解的改进算法
   • 原理：通过QR分解维持扰动向量的正交性，避免数值漂移和共线性问题，结合分数阶积分方法在每个时间步长对扰动向量进行演化和正交化。
   • 优势：适用于高维分数阶系统，计算稳定性高，可并行化处理。
   • 应用案例：
     • 分数阶Lorenz系统：通过改进算法计算Lyapunov指数谱，验证了分数阶阶次对混沌行为的影响。当阶次降低时，最大Lyapunov指数可能由正转负，系统从混沌状态转变为周期运动或稳定平衡点。
     • 4D分数阶Chen系统：揭示了超混沌行为（存在两个或多个正Lyapunov指数）与分数阶阶次的关联，阶次降低可能导致超混沌消失或吸引子结构变化。
2. 基于Wolf方法的推广
   • 原理：监测相空间中初始距离相近的轨道距离随时间的演化，定期重标化以避免数值溢出，通过局部扩张率估计Lyapunov指数。
   • 改进：结合分数阶系统的记忆效应，优化距离重标化策略，提升计算精度。
   • 应用案例：
     • 分数阶Duffing振荡器：分析分数阶阶次对混沌区域的影响，发现阶次降低可能抑制混沌或诱发混沌，取决于系统参数。

三、关键研究进展与成果

1. 分数阶Lorenz系统
   • 混沌行为与阶次关系：当分数阶阶次接近1时，系统行为与经典Lorenz系统相似；阶次降低时，混沌吸引子可能“收缩”或出现新的动力学行为（如周期运动）。
   • 数值验证：通过改进QR分解算法，计算结果与Wolf法对比，验证了算法有效性，并揭示了分数阶导数对混沌抑制的机制（记忆效应限制轨道发散）。
2. 4D分数阶Chen系统
   • 超混沌特性：当参数适当且阶次接近1时，系统存在两个正Lyapunov指数，表现为超混沌行为；阶次降低可能导致正指数减少，系统从超混沌转变为混沌或稳定状态。
   • 高维系统分析：4D结构为研究分数阶微积分对高维系统动力学的影响提供了平台，发现阶次变化可能改变相空间轨道维数和吸引子分形维数。
3. 分数阶Duffing振荡器
   • 阶次对混沌的调控：阶次接近1时，系统在参数混沌区域表现为混沌；阶次降低可能抑制混沌（轨道收敛）或诱发混沌（新吸引子出现），取决于阻尼比和策动力幅值。
   • 工程意义：为非线性振荡器的控制和同步设计提供了理论依据，如通过调节阶次实现混沌掩盖通信。

四、研究方法对比与推荐

五、未来研究方向

1. 算法优化：结合深度学习或自适应步长控制，进一步提升计算效率和精度。
2. 高维系统扩展：研究更高维分数阶系统的Lyapunov指数谱计算方法，如分数阶神经网络、电力系统模型。
3. 实验验证：通过物理电路或数值仿真平台（如Multisim、Simulink）验证理论结果，推动工程应用（如保密通信、混沌控制）。

六、结论

当前研究已建立适用于分数阶系统的Lyapunov指数谱计算方法，并通过分数阶Lorenz系统、4D Chen系统和Duffing振荡器验证了算法有效性。改进QR分解法和推广Wolf法为分数阶系统动力学分析提供了有力工具，未来需进一步优化算法并拓展应用场景，以揭示分数阶微积分在复杂系统建模中的潜力。

📚2 运行结果

2.1 测试1

2.2 测试2

2.3 测试3

部分代码：

% Output:
% x1,x2,x3,x4 - system responses;
% t - Time series of system responses;
% LE - Lyapunov Exponents;
% T - Time series of Lyapunov exponents;

%% system parameters and simulation conditions
clear;clc;close all;
C=0.3;
BETA=-0.1;
F=0.255;
W=1.2;
p=0.8;

h=1e-3;
h_norm=10*h;
N=h_norm/h;

tn=300-h;
t=0:h:tn;
n=length(t);
T=0:h_norm:tn;
%% define the order
q1=1;q2=1;q3=1-p;q4=1;

%% Fractional-order binomial coefficient
cp1=1; cp2=1; cp3=1; cp4=1;
for j=1:n
c1(j)=(1-(1+q1)/j)*cp1;
c2(j)=(1-(1+q2)/j)*cp2;
c3(j)=(1-(1+q3)/j)*cp3;
c4(j)=(1-(1+q4)/j)*cp4;
cp1=c1(j); cp2=c2(j); cp3=c3(j); cp4=c4(j);
end

%% initialization
x(1) = 0;
y(1) = 0;
z(1) = 0;
u(1) = 0;

🎉3参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。

🌈4 Matlab代码实现、数据、文章