Duffing系统共存周期解的稳定性与分岔演化

张锦涛，王 昕，吕小红，金 花

（兰州交通大学 机电工程学院，兰州 730070）

Duffing 系统是为描述机械系统中弹簧硬化效应而建立的一类典型非线性系统，很多工程中的振动问题可以使用该系统来研究。Duffing 系统虽然是比较简单的非线性系统，但是却包含丰富的非线性特性，如周期运动的跳变、周期解的共存、混沌激变等。目前，对于Duffing 系统的研究已有很多成果。文献[1]中提出了一种简单的周期参数切换方法，该方法可以找到广义Duffing系统中任何可以数值逼近的稳定极限环。文献[2]中研究了一类具有不连续非线性阻尼的Duffing方程周期解问题，证明其周期解只能是谐波解和次谐波解，然后利用庞加莱-波尔定理研究了该方程的解存在唯一性。文献[3]中研究了双势阱Duffing 方程混沌与周期吸引子的对称破缺激变。文献[4]中运用谐波平衡法计算系统的频率-振幅响应曲线，再通过预估-校正算法找到了系统随参数变化的跳跃区间，研究了其跳跃机理。文献[5]中采用数值方法计算双参数平面上的最大Lyapunov 指数，得到分岔曲线，分析了系统在双参数平面上的分岔及其全局特性。文献[6]中基于非线性振动系统的主要特征提出一种新的识别Duffing 非线性系统刚度的方法。文献[7]中应用谐波平衡法求解Duffing系统的振动方程，研究了常数激励与简谐激励联合作用下系统的骨架曲线及幅频响应特性。文献[8]中通过多尺度法得到系统的近似解析解，研究在多频激励下同时发生主共振和1/3次亚谐共振的动力学行为与稳定性。文献[9]中运用多尺度法求解一类含五次非线性恢复力的Duffing系统的幅频响应方程，分析系统的共振与分岔特性。以上文献多运用不同解析法对Duffing 方程的周期解进行计算，进而研究了系统中存在的多种非线性特性，但这些研究对于系统中共存周期解的产生原因及其分岔演化过程却鲜有提及。因此，系统中共存周期解的存在区域及其分岔演化规律没有完全被揭示。

本文将延续算法和打靶法相结合，延拓追踪系统中共存的周期解并应用Floquet 理论判断其稳定性及分岔。最后应用简单胞映射法计算系统共存周期解的吸引域，并结合相图和Poincaré 映射图研究周期解的吸引域随系统参数变化的侵蚀演变过程。

1 Duffing方程的延拓打靶法

式中，无量纲参数α为阻尼系数，β为刚度系数，ξ为非线性项系数，f为外激励幅值，ω为激励频率。为方便数值计算，令=v，则方程可降阶为：

为了研究系统参数变化时周期解的稳定性及其分岔演化规律，本文将延续算法和一般打靶法结合求解系统的周期解，并对其进行追踪分析[10]。

对于Duffing系统(2)，用X表示状态向量[x,v]T,用μ表示参数向量[α,β,ξ,ω]T,可简记为：

其中，求解系统的周期解需要构造系统的Poincaré映射，以此将周期解求解问题转化为满足下列条件的映射不动点曲线问题，即：

式中：P(μ,X)是X的Poincaré映射。

求解式(4)这样的非线性方程组，一般使用Newton-Raphson 方法进行迭代求解。延拓打靶法的基本思想是将参数空间[μ0,μ] 进行离散，当μ=μ0时求解式(4)可以得到X(μ0)，记为X0。将X0作为μ=μ1时的近似值X01，由式(4)进一步迭代可以得到X(μ1)，记为X1。以此类推得到一个近似的不动点X0k。此时可以用Newton-Raphson 方法进一步迭代求解，求得此时的不动点为Xk。由式(4)可得：

牛顿迭代格式如下：

针对部分老化管线结垢、腐蚀严重，注水生产时压力损失大、刺漏状况频发等问题，2017年完成17条腐蚀老化注水管线、1条注水干线的改造，将管线材质由复合管改为无缝钢管，改造总长度16.3 km，涉及17口注水井，对应注配间柱塞泵降压 0.5 MPa，节电 4.2×104kWh。

式中：i为迭代次数，实际计算中P(Xk)和DP(Xk)可由式(3)和式(7)以[X0k,I]为初始值积分一个Poincaré映射周期T得到。

为了提高算法的效率，对延续过程中的近似值先进行预估，再迭代校正。将式(4)转化为一般常微分方程的初值问题：

求解得到Xk时，运用Euler 法可以对μ=μk+1时的不动点进行预估得到X0k+1。

式中：Δμ为系统参数的变化量，∂G(μk,Xk)/∂X可由式(3)、式(5)、式(7)联合求解得到。∂G(μk,Xk)/∂μ可由式(3)和式(10)以[Xk,0]为初始值积分一个Poincaré 映射周期T得到。本文计算中选取激励频率ω为延拓参数，其余参数均为常数。

上述计算过程中，在求得不动点的同时得到了Poincaré 映射的Jacobi 矩阵DP(Xk)的特征值，根据Floquet 理论可判断周期解的稳定性。当所有Floquet 乘子都在复平面上的单位圆内时周期解稳定，存在一个Floquet乘子在单位圆外时周期解不稳定。随着参数的变化，当存在一个特征乘子从单位圆穿出时，周期解发生分岔而失稳，可从乘子穿出单位圆的方式识别其分岔类型。

2 基于延拓打靶法的分岔分析

Duffing 系统存在丰富的动力学行为。为研究周期解的分岔演化过程，取参数：α=0.83，β=1.0，ξ=0.5，f=2.0，数值计算系统的多初值分岔图如图1所示。

从图1中可以看出，在ω∈[0.8，0.965]时系统中只有稳定的周期一解P1a（灰线）。随着ω的增加，当ω=0.965 238 6 时，系统出现了新的稳定周期一解P1b（黑线），此时P1b的Floquet 乘子为(0.004 5，0.9997)，由Floquet 理论可知，P1b发生了鞍结分岔。在ω∈[0.965，1.032]内，稳定的周期一解P1a和P1b共存。随后，周期一解P1b经过叉式分岔产生两个稳定的反对称周期一解P1c和P1d。在ω∈[1.032，1.041]内，反对称的周期一解P1c和P1d以及周期一解P1a共存。当ω=1.040 942 6 时，周期一解P1a的Floquet乘子为(0.999 9，0.006 7)，由Floquet 理论可知，周期一解P1a发生鞍结分岔。当ω＞1.041 时，周期一解P1a消失，系统中只有反对称的周期一解P1c和P1d共存。P1c和P1d随着ω的增加经倍周期分岔序列进入混沌，此时，两个反对称的混沌吸引子共存。当ω=1.085 时，两个反对称混沌吸引子融合激变，生成一个全新的混沌吸引子。

图1 多初值分岔图

由图1可见稳定周期三解P3a的窗口。当ω减小时，P3a经叉式分岔产生两个反对称的周期三解P3b（灰线）和P3c（黑线），然后经倍周期分岔序列嵌入混沌。当ω增加时，P3a直接转迁为混沌。继续增加ω，当ω=1.833时，混沌吸引子再次激变为两部分，即产生两个反对称的混沌吸引子。随后混沌吸引子经倍周期倒分岔退出混沌，最终进入稳定的反对称周期一解的窗口。

为了深入研究系统周期解的稳定性与分岔，应用延拓打靶法追踪周期一解在ω∈[0.9，1.05]区间内的分岔演化。延拓计算之前，首先在系统的状态空间选择一个考察区域Ω=(x，v)=(-3≤x≤3，-3≤v≤3)，然后将状态平面Ω划分为100×100=10 000 个小网格。以每个小网格的中心点为初始值应用打靶法计算系统的周期解，则共有10 000 个初始值需要考察。当ω=1 时，求得系统有3 个共存的周期解（2 个稳定的周期一解(0.031 239 28，1.746 523 36)和(-1.028 217 48，0.657 523 80)，1个不稳定的周期一解(-1.066 474 56，1.417 145 48)）。分别以每个周期解为初值进行延拓计算，结果如图2所示。图中实线表示稳定的周期解，“+”标记线表示不稳定的周期解。

图2 周期一解的延拓追踪

由图2可知，在ω∈[0.965，1.032]内，系统存在两个稳定的周期一解P1a和P1b，以及一个不稳定的周期一解UP1a。在ω∈[1.032，1.041]内，3 个稳定的周期一解P1a、P1c和P1d，以及两个不稳定的周期一解UP1a和UP1b共存。在ω∈[1.041，1.05]内，两个稳定的周期一解P1c和P1d，以及一个不稳定的周期一解UP1b共存。当ω=1.041时，稳定的周期一解P1a与不稳定的周期一解UP1a相撞后消失，由前面的分析可知，此时系统发生鞍结分岔。当ω=0.965时，稳定的周期一解P1b与不稳定的周期一解UP1a经鞍结分岔消失。由此可见，鞍结分岔是Duffing系统中共存吸引子消失的重要原因。当ω=1.032 时，稳定的周期一解P1b失稳产生不稳定的周期一解UP1b，同时通过叉式分岔产生两个稳定的反对称周期一解P1c和P1d。

以ω=1.41时的8个共存周期三解为初值应用延拓打靶法延拓追踪每个周期三解在ω∈[1.41，1.70]内的分岔演化如图3所示。图中稳定的周期解用实线表示，不稳定的周期解用不同符号的虚点区分。由图3可知，在ω∈[1.444，1.632]内，一个稳定的周期三解P3i与5个不稳定的周期三解P3d至P3h共存。图3(b)为图3(a)黑框中不稳定周期三解的局部放大。当ω减小至1.444 时，周期三解P3i失稳转迁为不稳定的周期三解P3c，同时通过叉式分岔产生两个稳定的反对称周期三解P3a和P3b，在此分岔过程中，周期三解P3d至P3h的模式及稳定性保持不变。因此，在ω∈[1.41,1.444]内，系统存在两个稳定的反对称周期三解P3a和P3b，以及6个不稳定的周期三解P3c至P3h。图4(a)至图4(h)依次给出了ω=1.41 时周期三解P3a至P3h的相图和Poincaré 映射图。当ω增加至1.631 978 7时，稳定的周期三解P3i与不稳定的周期三解P3d碰撞之后消失，此时P3i的Floquet 乘子为(0.999 6，0.000 1)，结合图1并由Floquet 理论可知，周期三解P3i通过鞍结分岔直接进入混沌。因此，在ω∈[1.632，1.700]内，4个不稳定的周期三解P3e至P3h与混沌吸引子共存。在周期三解的窗口中，由于稳定的周期三解和多个不稳定的周期三解共存，导致系统的稳定性降低。

图3 周期三解的延拓追踪和局部放大

图4 ω=1.41时共存周期三解的相图和Poincaré映射图

图4 ω=1.41时共存周期三解的相图和Poincaré映射图

由以上分析可知，Duffing系统在分岔图的周期解窗口中，稳定和不稳定周期解的共存现象普遍存在，而不稳定周期解经过较长时间后会趋于稳定的周期解。当稳定周期解共存时，系统的响应与初值条件密切相关。采用周期解的吸引域能够更为清晰地描述系统周期的初值响应和全局特性。因此，下文将根据简单胞映射法计算稳定周期解的吸引域，进一步揭示系统随参数变化时周期解吸引域的演化规律以及系统的全局分岔特性。

3 共存吸引子的吸引域演化

当多个吸引子共存时，系统的动力学行为受初始条件的影响较大。选用第2 节所用参数，对系统的全局动力学特性进行研究。在系统的相空间中选取一个考察区域H={-4.0≤x≤4.0，-4.0≤v≤4.0}，将其划分为800×800 个胞。以每个胞的中心点为初值，应用简单胞映射法求解图2中稳定周期一解的吸引域，揭示激励频率ω增大时系统周期解的吸引域演化与分岔过程，如图5所示。

由第2 节中的分析可知，在ω∈[0.9，0.965]时系统只存在一个稳定的周期一解P1a。图5(a)为ω=0.950 时的吸引域，其中白色“×”号表示周期一吸引子P1a，黑色区域为其吸引域，此时系统在考察区域H 内只存在一个周期解。对应的相图和Poincaré 映射图如图6(a)所示。

随着ω增大，当ω=0.965 时系统经鞍结分岔产生新的稳定周期一解P1b，并与P1a共存。图5(b)给出了ω=0.968时共存吸引子的吸引域，其中白色“×”号表示周期一吸引子P1a（94.09%），白色“+”号表示周期一吸引子P1b，灰色区域为P1b的吸引域（5.91%）。对比图5(a)，可见在H 内P1a原本完整的黑色吸引域中出现了P1b的灰色吸引域，即P1a在部分初值范围内发生周期跳跃转迁为P1b，P1a的吸引域被P1b的吸引域侵蚀。当ω=1.0时共存吸引子P1a和P1b的吸引域如图5(c)所示，可见随着ω增大，系统中灰色吸引域的面积不断增大，P1a的吸引域被进一步侵蚀。此时P1b的稳定性增强，但是其吸引域面积（15.79%）远小于P1a的吸引域面积（84.21%），P1a的稳定性仍然最强。P1a和P1b对应的相图和Poincaré映射图如图6(b)所示，两运动轨迹的振动幅值不同，P1b的振幅相对较小。进一步增加ω至1.025时，P1a（85.71 %）和P1b（14.29 %）的吸引域如图5(d)所示。由图可见，P1a在H内的右侧被P1b部分侵蚀，而P1b集中分布的灰色吸引域被P1a侵蚀，即原来的P1b在部分初值范围内转迁为P1a，两个吸引域嵌套在一起，其初始条件的微小扰动有可能使系统出现不同的周期响应，P1b的稳定性明显降低。

继续增大ω，当ω=1.032时，P1b经叉式分岔产生两个反对称的周期一吸引子P1c和P1d。当ω＞1.032时，系统中3个吸引子P1a、P1c和P1d共存。图5(e)给出了ω=1.035时不同吸引子的吸引域，其中“×”号表示P1a，黑色为其吸引域（65.22%），“◇”号表示P1c，灰色为其吸引域（19.92%），“△”号表示P1d，白色为其吸引域（14.86%）。图中两个反对称吸引子P1c和P1d的吸引域互相嵌套并缠绕在一起，两个吸引子的稳定性较弱，而P1a的吸引域面积较大，稳定性较强。3个共存吸引子的相图和Poincaré映射图如图6(c)所示。随着ω增大，P1a的吸引域（56.76%）被缠绕在一起的P1c（21.27%）和P1d（21.97%）的吸引域分别从下往上和从右往左逐渐侵蚀，P1a的稳定性逐渐降低，如图5(f)所示。

当ω增加至1.041时，P1a发生鞍结分岔而消失，系统中只剩下P1c和P1d共存。ω=1.042 时，两个共存吸引子的吸引域如图5(g)所示。由图可见，系统在H 内只存在P1c（48.99%）和P1d（51.01%）的吸引域，原本P1a的吸引域被完全侵蚀。图5(h)给出了ω=1.05 时两个反对称吸引子P1c（50.81 %）和P1d（49.19%）的吸引域。对比图5(g)可见，随着ω增大，P1c的吸引域和P1d的吸引域面积逐渐接近，并且缠绕在一起，具有相似的分形结构。此时系统对初值极度敏感，初值的微小变化会使系统运动到不同的周期轨道上。P1c和P1d对应的相图和Poincaré映射图如图6(d)所示。

图5 随激励频率ω增大吸引域的演化与岔过程

图6 相图和Poincaré映射图

由以上分析可知，系统中出现共存吸引子时，随着ω增大，新吸引子的稳定性会逐渐增强，如图5(b)和图5(c)所示，P1b的吸引域从中心向外扩展，面积逐渐增大；但是在分岔点附近吸引子的稳定性会减弱，如图5(d)和图5(f)所示，P1b和P1a在发生分岔前原本完整的吸引域被与其共存的其余吸引子的吸引域所侵蚀。鞍结分岔会引起系统新的周期解的出现或消失，从而导致考察区域内周期解类型的突变及吸引域拓扑结构的改变，其稳定性也随之发生变化，如图5(a)和图5(b)以及图5(f)和图5(g)所示。当系统中由于叉式分岔出现、反对称吸引子共存时，反对称吸引子的吸引域拓扑结构相似，互相缠绕，导致系统的稳定性降低。吸引域的演化和侵蚀对系统周期解的稳定性有重要的影响。

4 结语

本文应用打靶法和Floquet 理论计算得到Duffing系统稳定和不稳定的周期解，并对其进行延拓追踪，分析了共存周期解的分岔演化过程及其全局动力学特性。得到了以下结论：

(1)Duffing 系统中多个稳定及不稳定周期解的共存现象较为普遍，并具有一定反对称性。当反对称吸引子通过倍周期分岔进入或退出混沌时，会发生混沌吸引子的内部激变。

(2)稳定周期解与不稳定周期解碰撞而发生鞍结分岔是共存吸引子出现和消失的重要原因，并且鞍结分岔会导致吸引子的吸引域拓扑结构发生改变。当系统由于叉式分岔产生、两个反对称吸引子共存时，其吸引域拓扑结构相似，彼此嵌套缠绕在一起，会导致系统的稳定性降低。Duffing系统的周期吸引子主要通过倍周期分岔或鞍结分岔进入混沌。

(3)当系统中出现共存吸引子时，稳定周期解受初始条件的影响较大。随着系统参数的变化，新出现吸引子的稳定性会逐渐增强，但在分岔点附近，即将发生分岔的吸引子稳定性会降低。