基于Melnikov方法的分数阶Duffing振子混沌阈值解析研究

秦 浩，温少芳，申永军，邢海军，王 军

(1.石家庄铁道大学 交通运输学院 省部共建交通工程结构力学行为与系统安全国家重点实验室，石家庄 050043；2.石家庄铁道大学 机械工程学院，石家庄 050043)

分数阶微积分已经存在了300多年，但是目前还是一个比较新鲜的概念。早在整数阶微积分创立初期，法国学家Hospital和德国数学家Leibniz就提出了分数阶的概念[1-3]，但是由于缺乏实际应用背景的支撑等诸多方面原因使得它长期以来没有得到研究和发展。随着自然科学、社会科学的发展和工程的实际需要，各种复杂的系统和计算机技术引入到工程领域，分数阶微积分的定义、特性和计算等才得以在工程领域应用。目前，分数阶微积分在工程领域中的应用主要分为两类：一类是分数阶微积分的控制，能够提高控制的效果和系统的鲁棒性；另一类是利用分数阶微积分建模，很多工程材料介于理想固体与理想流体之间，具有记忆特性，用分数阶微积分来模拟更加准确，更能反应真实的本构关系[4-8]。

早期学者在解释生活中的经典现象时提出了许多假说，其中大部分是基于整数阶模型模拟自然现象，并且模拟后的结果能够满足当时的需求，因此人们忽视了分数阶模型。随着研究的深入，研究者们发现对于一些具有黏弹性特性的材料单纯用整数阶的模型无法满足精度的需要，工程系统中的非光滑性、不连续性、参数激励、间隙等，使得系统响应变得更加复杂，并且会降低系统的整体性能，甚至会导致系统出现分岔和混沌的动力学现象[9-14]。因此人们将整数阶模型进一步优化为分数阶模型，用分数阶微分方程来模拟非经典现象，能够更精确地反应材料的本构关系，更好地满足实际的需求。在控制方面，分数阶相较于整数阶的优点在于能够控制或改变系统的闭环控制特性，增强控制效果，提高系统的鲁棒性。所以大量的学者开始研究分数阶微分方程的典型力学特性和分数阶微积分对整个动力系统力学特性的影响[15]。

随着分数阶微积分的发展，分数阶微积分方程的求解成为学者们研究的一个热点方向。分数阶微积分的计算较为复杂，大部分的文献中对分数阶的处理采用数值计算的方法，除此之外还有定性和解析的方法。数值计算方法是处理分数阶微积分的一个主流分支，应用较为广泛，如Diethelm等[16]讨论了一种可以应用于分数阶微分方程的求数值解的Adams型预测-矫正法；Firdous等[17]提出了一种基于哈尔小波求解分数阶微分方程的新的运算矩阵方法用于求解线性和非线性分数阶微分方程的数值结果；Javidi等[18]提出了一种基于同伦微扰法和拉普拉斯变换求解分数阶偏微分方程的数值方法。定性研究则侧重分析方程解的数目和稳定性的变化，比如Trigeassou等[19]利用李雅普诺夫方法通过频率分布分数阶微积分模型研究了分数阶微分方程的稳定性；Qi等[20]提出了分数阶的阶次大于1小于2时Caputo分数阶导数的一个新定理和Caputo分数阶系统的两个新定理，方便了使用李雅普诺夫方法证明分数阶非线性时变系统的稳定性。解析研究是通过求得近似解析解对方程进行定量分析，比如Shen等[21-23]利用平均法分别研究了分数阶Duffing振子和分数阶van der Pol振子的近似解析解，提出了等效线性刚度和等效线性阻尼的概念，分析了分数阶微分的参数对不同振子动力学特性的影响；姜源等[24]采用平均法研究了含分数阶微分项的Duffing振子的超谐与亚谐联合共振，得到了一阶近似解析解，分析了分数阶微分的参数对系统动力学特性的影响，并且提出了超、亚谐联合共振时等效线性阻尼和等效线性刚度的概念；顾晓辉等[25]利用多尺度法得到了两个谐波激励作用下分数阶Duffing振子的一类组合共振的一阶近似解析解，并且研究分析了定常解的稳定性；Tabejieu等[26]通过平均法研究了分数阶阶次对非线性Rayleigh梁振幅幅值的影响，证明了随着分数阶阶次的增加，梁的共振振幅减小，同时利用Melnikov方法分析了混沌产生的必要条件。

关于Duffing振子产生混沌必要条件的问题，很多学者利用Melnikov方法对其进行了研究。如：Shen等[27]基于Melnikov方法分析了带有位移延迟反馈和速度延迟反馈的整数阶Duffing振子产生混沌的必要条件，分析了位移延迟反馈和速度延迟反馈分别影响混沌产生必要条件的规律。Sun等[28]利用Melnikov方法研究了具有位移延迟反馈和速度延迟反馈的整数阶Duffing系统在谐波激励下的混沌行为，导出了同宿分岔产生混沌的必要条件的解析解。Wen等[29]用Melnikov方法研究了强迫激励下的整数阶Duffing振子在位移延迟反馈和速度延迟反馈控制下的异宿分岔和混沌现象。张思进等[30]研究了具有立方非线性项和外部激励项的二自由度非线性碰振系统的动力学特性，分别运用改进后的局部亚谐Melnikov方法和数值方法得到了二自由度碰振系统稳定周期运动的存在条件。Abtahi[31]采用Melnikov方法研究了陀螺卫星复杂自旋轨道动力学中的混沌现象。多数文献主要针对整数阶Duffing振子产生混沌的问题进行研究，关于分数阶Duffing振子产生混沌必要条件多采用数值解，解析结果还未见相关文献。下面重点关注应用Melnikov方法研究含有分数阶微分项的Duffing振子的分岔和混沌现象的问题。Melnikov方法是一种预测系统混沌的一阶近似方法，本文基于Melnikov方法提出一种分数阶Duffing振子产生混沌必要条件的新的研究思路。首先将分数阶微分项进行了处理，利用等效阻尼和等效刚度的概念，将分数阶微分项等效成三角函数与指数函数的形式。接着利用Melnikov方法研究分数阶Duffing振子的同宿轨分岔与混沌的产生条件。最后将解析结果与数值结果进行了比较，验证了解析结果的准确度，并分析了分数阶的阶次和系数对分数阶Duffing振子产生混沌的必要条件的影响。

1 分数阶Duffing振子产生混沌的必要条件

研究如下的含分数阶微分项的Duffing振子

(1)

(2)

式中，Γ(n)为Gamma函数，Γ(y+1)=yΓ(y)。

设式(1)的解满足

x(t)=acos(ωt+θ)

(3a)

(3b)

利用Caputo定义求解分数阶微分项，将式(3)代入式(2)，得到

(4)

引入s=t-u，则式(4)变为

(5)

引入两个基本公式[32]

(6a)

(6b)

代入式(5)，得到

(7)

代入方程原参数，得到

(8)

本文主要研究式(1)经过长时间振动稳定以后的动力学现象，所以在后面的研究中将高阶项略掉，只取一阶近似项，即

(9)

式(9)得到的结果与很多文献中得到的结论相同，分数阶微分项不仅起着阻尼的作用，同时也有刚度的作用，可以等价为等效线性阻尼和等效线性刚度。下面利用等效线性刚度与等效线性阻尼的概念，将式(9)代入式(1)变为

(10)

整理式(10)可得

(11)

引入c=εc1,λ=ελ1,f=εf1。

其中，ε为小参数，则式(11)变为

(12)

当ε=0时，则系统为未扰动系统

(13)

在式(13)中，存在一个同宿轨道满足

(14)

(15)

对式(14)进行积分计算，得到

(16)

计算式(16)的定积分，整理得到

(17)

将式(12)转化为

(18a)

其中，

(18b)

(18c)

(18d)

同宿轨道的位移为

(19a)

(19b)

建立Melnikov函数[33]

(20)

将式(20)中的定积分分为两部分M1(t0)和M2(t0)，分别经过定积分计算，得到

(21a)

根据函数的奇偶性质，式(21a)变为

(21b)

(22)

联立式(20)～式(22)，可以得到small意义下混沌的必要条件为

(23a)

代入原系统参数，式(23a)可变为

(23b)

当分数阶微分项产生的等效刚度小于系统的固有刚度时，将式(23b)中的绝对值去掉，得到分数阶Duffing系统产生混沌的必要条件为

(24)

2 数值仿真

为了验证本文解析结果的正确性，下面将数值结果与本文解析结果进行对比。数值仿真利用Diethelm等研究中的预测矫正法对分数阶Duffing系统式(1)进行计算，将计算结果中不稳定的响应去掉，选取稳态响应的部分，分析系统的稳态响应情况。

本文中选取的基本参数c=0.2，k=1，α=1，ω=1.1，λ=0.1，分数阶微分项的阶次以p=0.5为例进行仿真验证。系统初值首先取为“1，0，0，0”，计算时间取为10 000 s，舍去前5 000 s，确保数值的稳定性。并使振幅f由0～0.5变化，选择步长为0.005，得到系统式(1)产生分岔到混沌的过程示于图1(a)中。从图1(a)中可以看出系统产生混沌的临界外激励振幅fmin约为0.23，放大f=0.23周围的分岔图见图1(b)，横轴为外激励振幅f从0.20～0.26，步长选取0.001。

从图1(b)中可以看出混沌产生的临界振幅fmin约为0.234。为了证明f=0.234为产生混沌的临界点，下面f分别取为0.21，0.22，0.227，0.233，0.234时，分析其系统响应的特性如下：

(1) 当f=0.21时，系统的相图和时间历程图如图2(a)所示，此时的系统运动状态为单周期运动。

(2) 当f=0.22时，系统的相图和时间历程图如图2(b)所示，此时的系统运动状态为周期2的运动。

(3) 当f=0.227时，系统的相图和时间历程图如图2(c)所示，此时的系统运动状态为周期4的运动。

(4) 当f=0.233时，系统的相图和时间历程图如图2(d)所示，此时的系统运动状态为拟周期的运动。

图2 相图与时间历程图Fig.2 Phase portrait and time history

(5) 当f=0.234时，Duffing系统的相图和时间历程图如图3所示，从图中可以看出系统此时处于混沌状态。

由图1～图3的仿真结果，可以看出：当p=0.5时，外激励振幅f=0.234是利用数值方法得到的分数阶Duffing振子产生混沌的最小值。再根据本文解得的产生混沌必要条件的解析式(24)计算，得到系统产生混沌的必要条件的解析结果为fmin=0.22。对比数值解与解析解，差值为0.014。Melnikov方法本身是一种求解产生混沌必要条件的一阶近似方法，所以数值解与解析解必然会存在着误差。为了判断差值是否在合理范围内，将分数阶微分项设为零，其他参数取值同上。使用龙格库塔法求解此整数阶Duffing系统产生混沌的fmin数值解为0.192，而利用Melnikov方法计算的整数阶Duffing振子产生混沌的必要条件fmin为0.173，整数阶情况下数值解与解析解的差值为0.019。两个差值属于同一量级，说明这种误差主要源自Melnikov方法本身。而本文中提出的利用等效阻尼和等效刚度将分数阶系统等效替换成整数阶系统进行分析，精确度较高，方法可行。

图1 分岔图Fig.1 Bifurcation diagram

图3 相图与时间历程图Fig.3 Phase portrait and time history

在数值仿真过程中，还发现改变系统的初值时，出现了另外一条通过倍周期分岔通往混沌的路径，如：当初值取为“-1，0，0，0”时，其他系统参数的取值同上，得到的系统分岔图如图4所示。当f分别为0.21，0.22，0.227，0.233和0.234时，系统的相图和时间历程图如下：

(1) 当f=0.21时，系统的相图和时间历程图如图5(a)所示，此时的系统运动状态为单周期运动。

(2) 当f=0.22时，系统的相图和时间历程图如图5(b)所示，此时的系统运动状态为周期2的运动。

(3) 当f=0.227时，系统的相图和时间历程图如图5(c)所示，此时的系统运动状态为周期4的运动。

(4) 当f=0.233时，系统的相图和时间历程图如图5(d)所示，此时的系统运动状态为拟周期的运动。

(5) 当f=0.234时，分数阶Duffing系统的相图和时间历程图如图6所示,与图3(a)和图3(b)运动状态相同，均为混沌状态。

通过图1和图4、图2和图5、图3和图6的两两对比分析可知：当f取值小于0.234时，选取上述不同初值时，分数阶Duffing振子的运动特性都是相同的，相图方向相反，分别为图3和图6中的混沌相图的左右半支。但是不同初值时，产生混沌的临界位置并没有改变，fmin仍为0.234。所以，可以证明改变系统的初值不会改变系统产生混沌的必要条件,但是取不同初值分数阶Duffing振子会通过不同的倍周期分岔路径到达混沌。主要原因是该分数阶Duffing振子具有双稳态特性，在文中参数条件下，外激励f取值小于图4所示的分岔点f1处时，系统处于单周期运动，存在两个稳态解。取不同初值时，将进入不同的吸引域从而形成不同的稳态解。从两个稳态解出发，随着外激励参数的变化都能通过倍周期分岔到达混沌的状态。

图4 p=0.5的分岔图Fig.4 Bifurcation diagram for p=0.5

图5 相图与时间历程图Fig.5 Phase portrait and time history

图6 相图与时间历程图Fig.6 Phase portrait and time history

3 分数阶参数的影响分析

下面主要研究分数阶微分项的阶次和系数对分数阶Duffing振子产生分岔和混沌的必要条件的影响。

首先分析分数阶微分项的阶次对系统产生混沌的必要条件的影响。选取的基本参数为c=0.2，k=1，α=1，ω=1.1，λ=0.1，分数阶阶次p在[0 1]取值，利用本文解得的解析结果，即式(24)，得到系统产生混沌必要条件的激励振幅fmin与分数阶阶次p的关系曲线，示于图7中用实线表示。从图7可以看出，当分数阶阶次p从0增大到1的过程中，系统产生混沌的必要条件fmin是逐渐增大的。为了验证本文的解析结果的精确度，在不同p值时，利用数值仿真的方法也进行了验证，图中圆圈分别代表不同分数阶阶次下的数值仿真结果，从两者之间的对比可以看出有一定的差值，本文得到的解析值在定量上小于数值结果，但整体趋势相同，定性分析是一致的。可见分数阶阶次的增加能抑制系统的分岔和混沌的产生。

图7 数值解和解析解比较Fig.7 The comparisons between the numerical and analytical solutions

接着分析分数阶微分项的系数对系统产生混沌的必要条件的影响。选取的基本参数为c=0.2，k=1，α=1，ω=1.1，分数阶阶次p在[0 1]取值，如图8(a)所示。经过对比发现改变分数阶项系数λ变大时，系统产生混沌的必要条件要求fmin也要不断增大，可见分数阶微分项的系数的增大同样也会抑制分数阶Duffing振子产生分岔和混沌的现象。

选取分数阶阶次p为0.5时，系统产生混沌的必要条件fmin随着分数阶系数改变的变化规律见图8(b)所示，从图中可以看出，随着分数阶微分项系数λ的增加，fmin也会逐渐增大，当分数阶微分项系数λ远远小于线性刚度k时，解析解和数值解之间的差别较小。当分数阶微分项系数λ继续增大，分数阶微分项的等效线性刚度接近于系统的线性刚度k时，fmin出现增大较快的现象，此时解析解和数值解之间的差值会加大。这是因为利用等效线性阻尼和等效线性刚度来替代分数阶微分项的方法省略了式(8)的高阶项，当分数阶微分项的等效线性刚度与线性刚度k越接近时，省略的两项高阶项所占的比值就会越大，所以此时解析解和数值解差值会加大。因此本文的解析方法要求分数阶微分项的系数要远小于系统的线性刚度k。

图8 分数阶系数的影响Fig.8 The influence of fractional order coefficients

4 结 论

本文利用Melnikov方法研究了含有分数阶微分项的Duffing振子产生混沌的必要条件，将分数阶Duffing振子中分数阶微分项替换成等效线性阻尼与等效线性刚度，用Melnikov方法分析等效后的整数阶系统的同宿轨分岔和产生混沌的必要条件，并得到了其解析结果，并通过数值迭代的方法验证了其精确度。通过对比不同阶次情况下数值解与解析解得到的产生混沌的阈值，发现解析解的趋势与数值迭代模拟一致，临界激励振幅fmin会随着分数阶阶次的增加而增大，定性结果相同，Melnikov方法分析系统混沌产生的必要条件是一阶近似结果，因此数值结果与解析解之间的定量差异是可以接受的。最后分析了分数阶微分项的系数对系统产生混沌的必要条件的影响。在数值模拟过程中,还发现在分数阶Duffing振子中，当选取不同的初值时，存在两条通过倍周期分岔通往混沌的路径，并通过分析其动力学响应证实了这一现象。本文提出的将分数阶微分项替换为三角函数和指数函数形式的等效线性阻尼和线性刚度的方法可以推广应用于其他类似的分数阶系统中，为分数阶系统的动力学研究提供了新的研究思路。