浙江师范大学数学与计算机科学学院，浙江 金华

收稿日期：2020年12月28日；录用日期：2021年1月20日；发布日期：2021年1月29日

摘要

本文主要探究一类有周期激励的Duffing系统，根据Melnikov方法、欧拉离散、四阶Runge-Kutta算法理论及Matlab程序等进行数值模拟。从理论上分析不同参数下系统的动力学现象，如分岔、周期运动、混沌运动等，作出时间序列图、相图、分岔图及最大Lyapunov指数图，观测在振幅、频率及阻尼项变化时Duffing系统产生的复杂动态行为。数值模拟发现，在考虑激励振幅和频率不变时，随阻尼项的增大，系统混沌行为会发生显著变化，混沌状态将逐渐消失形成周期运动，而在考虑振幅影响时，系统混沌状态则无明显改变。

关键词

Duffing系统，阻尼系数，混沌现象，四阶Runge-Kutta算法

Numerical Simulation of a Class of Duffing Systems with Periodic Excitation

Liang Li, Youhua Qian^*

College of Mathematics and Computer Science, Zhejiang Normal University, Jinhua Zhejiang

Received: Dec. 28^th, 2020; accepted: Jan. 20^th, 2021; published: Jan. 29^th, 2021

ABSTRACT

This paper mainly explores a Duffing system with periodic excitation. Based on Melnikov method, Euler discretization, fourth-order Runge-Kutta algorithm theory and Matlab program, it theoretically analyzes the dynamic phenomena of the system under different parameters, such as bifurcation, periodic motion and chaos, makes time series diagram, phase diagram, bifurcation diagram and maximum Lyapunov index diagram, and observes Duffing system when amplitude, frequency and damping terms change. Numerical simulation shows that when the excitation amplitude and frequency are constant, with the increase of damping term, the chaotic behavior of the system will change significantly, and the chaotic state will gradually disappear to form periodic motion, while when the influence of amplitude is considered, the chaotic state of the system will not change much.

Keywords:Duffing System, Damping Coefficient, Chaotic Phenomenon, Fourth-Order Runge-Kutta Algorithm

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

随着科学技术的发展，非线性问题 [1] 出现在越来越多的学科之中。对于非线性振动系统的研究，大多建立于各种物理学模型，如Lorentz模型、悬臂梁模型、单摆的有限振幅振动模型、双自由度系统振动模型等。在非线性振动 [2] 理论研究中，Duffing方程 [3] 占有重要地位，对于含有周期激励的系统而言，由于系统的固有频率和外加周期激励项的频率相互作用，使得系统蕴含丰富的动态响应行为，如混沌 [4]、倍周期分岔、簇发振荡 [5] 等现象。特别是具有非线性恢复力的系统受到简谐激励 [6] [7] 时，受迫振动 [8] 存在跳跃现象，即激励频率或阻尼参数在通过某些特定值时，系统将产生跳跃突变或混沌。

在混沌控制 [9] 范畴下，对于目前确定性的混沌系统 [10]，还较少考虑各参数激励的影响，即在阻尼和外激励作用下，如何设置有效参数范围控制系统的混沌行为。对于含参数的非线性系统而言，其稳定性不仅取决于控制系统的固有结构和参数，而且与系统的初始条件以及外加输入有关，因此需要参数控制来对系统可能产生的动态行为进行模拟分析。针对非线性Duffing方程，通过数值计算 [11] 和数值模拟来刻画动力学行为，并依据已有的物理、力学振动模型 [12] 来反映实际问题中的混沌现象。因此，在非线性动力学系统中，有效处理参数项对系统具有重要意义。

本文应用Melnikov方法 [13]、欧拉离散、四阶Runge-Kutta算法 [14] 理论、混沌理论和分岔理论 [15] [16]，动态模拟一类具有周期激励Duffing方程随系统参数变化时的复杂动力学行为，以及各参数对系统混沌行为的影响，这些研究将有助于我们进行下一阶段的理论分析和模拟实验。

2. Duffing系统的混沌机理及离散化

这部分主要分析这一类具有周期激励Duffing系统，通过哈密顿函数、Melnikov法计算异宿轨线及Melnikov阈值 [17]，来确定系统存在的混沌现象，并用欧拉方法对系统进行离散，动态模拟不动点分岔。

2.1. 考虑非线性Duffing系统

$x^{″}+2\mu \alpha x^{\prime }+mx-n{x}^3=\mu f\cos \Omega t$ (1)

其中，阻尼系数 $\mu >0\left(0<\mu <1\right)$，外激励振幅f，频率 $\Omega$，非线性刚度系数为n，且 $m=\delta \cos \alpha t$，非线性恢复力为 $-mx+n{x}^3$，那么在振幅值 $\mu f$ 和激励频率 $\Omega$ 周期外力作用下系统会表现出周期或混沌行为，且在外激励频率 $\Omega \ll \alpha$ 时，系统存在多尺度的复杂动力学行为。

现将上式转化为二维非线性微分方程如下：

$\begin{array}{l}{x}^{\prime }=y\\ y^{\prime }=-2\mu \alpha y-mx+n{x}^3+\mu f\cos \Omega t\end{array}$ (2)

令 $\mu =0$ 时，考虑 $m=n=1$，上式的哈密顿系统对应的函数如下：

$H\left(x,y\right)=-\frac{1}{4}{x}^4+\frac{1}{2}{x}^2+\frac{1}{2}{y}^2$ (3)

该系统存在三个平衡点为两个鞍点 $\left(\pm 1,0\right)$，一个中心点 $\left(0,0\right)$，进而得到两条异宿轨如下：

$x_{\pm }^0\left(t\right)=\pm \tanh \left(0.5\sqrt{2}t\right);y_{\pm }^0\left(t\right)=\pm 0.5\sqrt{2}{\text{sech}}^2\left(0.5\sqrt{2}t\right)$ (4)

现用Melnikov方法分析异宿轨破裂产生混沌现象的条件，因此得出Melnikov函数为：

$\begin{array}{l}{M}_{\pm }\left(t_0\right)={\displaystyle {\int }_{-\infty }^{+\infty }{y}_{\pm }^0\left(t\right)\cdot \left[-2\alpha y_{\pm }^0\left(t\right)+f\cos \Omega \left(t-t_0\right)\right]\text{d}t}\\ =-2\alpha {\displaystyle {\int }_{-\infty }^{+\infty }{\left(y_{\pm }^0\left(t\right)\right)}^2\text{d}t}+f{\displaystyle {\int }_{-\infty }^{+\infty }{y}_{\pm }^0\left(t\right)}\cos \Omega \left(t-t_0\right)\\ =A\alpha +f{\displaystyle {\int }_{-\infty }^{+\infty }{y}_{\pm }^0\left(t\right)}\left(\cos \Omega t\cos \Omega t_0+\sin \Omega t\sin \Omega t_0\right)\\ =A\alpha +B\cdot f\cos \Omega t_0\end{array}$ (5)

其中 $A=\frac{2\sqrt{2}}{3},B=\sqrt{2}\pi \Omega \text{csch}\left(0.5\sqrt{2}\pi \Omega \right)$ 此时，当 $\frac{f}{\alpha }>\frac{A}{B}$ 时，即：

当 $\frac{f}{\alpha }>\frac{2\sqrt{2}}{3}\cdot \frac{1}{\sqrt{2}\pi \Omega \text{csch}\left(\text{0}\text{.5}\sqrt{2}\pi \Omega \right)}$，存在 $t_0^{\pm }$，使得 $M_{\pm }\left(t_0^{\pm }\right)=0$，又因为 $M_{\pm }\left(t_0^{\pm }\right)=-\Omega fB\sin \left(\Omega t_0^{\pm }\right)\ne 0$，所以 $M_{\pm }\left(t_0\right)$ 存在单重简单零点，则会产生Smale马蹄意义下的混沌现象。其中 $\frac{A}{B}=\frac{2}{3\pi \Omega }\sinh \left(0.5\sqrt{2}\pi \Omega \right)$，因此，当外激励振幅f与 $\alpha$ 之比超过Melnikov阈值 $\frac{A}{B}$ 时，Duffing系统会产生混沌现象。

此时不妨考虑 $\Omega =0.1$，则混沌阈值 $\frac{A}{B}=0.467$，即 $\frac{f}{\alpha }>0.467$ 时产生混沌现象，则选取 $f=0.35,\alpha =0.7$，有 $\frac{f}{\alpha }=0.5>\frac{A}{B}$，并取 $\mu =0.1$，对系统进行数值模拟，如图1表示该系统混沌状态下的相图和时间序列图。

图1. 当 $\mu =0.1,f=0.35,\alpha =0.7,\Omega =0.1$ 时，系统相图和时间序列图

2.2. 对非线性Duffing系统离散化

对于 $\mu \ne 0$，考虑方程(2)设 $p=\mu f\cos \Omega t$，利用欧拉方法，得出二维离散系统为：

$\begin{array}{l}{x}_{n+1}=x_n+\Delta t{y}_n\\ y_{n+1}=-m\Delta t{x}_n+\left(1-2\mu \alpha \Delta t\right)y_n+n{\left(\Delta t\right)}^3{x}_n^3+\Delta tp\end{array}$ (6)

设 $k_n=x_n+\Delta t{y}_n$，方程(6)变为：

$\begin{array}{l}{x}_{n+1}=k_n\\ k_{n+1}=2\left(1-\mu \alpha \Delta t\right)k_n+\left(2\mu \alpha \Delta t-m\Delta t^2-1\right)x_n+n\Delta t^4{x}_n^3+\Delta t^{2}p\end{array}$ (7)

接着，对方程(7)作新的尺度变换为：

$\begin{array}{l}{x}_{n+1}=y_n\\ y_{n+1}={\delta }_1{x}_n+{\delta }_2{x}_n^3+{\delta }_3{y}_n+p\end{array}$ (8)

其中， ${\delta }_1=2\mu \alpha \Delta t-m\Delta t^2-1,{\delta }_2=n\Delta t^2,{\delta }_3=2\left(1-\mu \alpha \Delta t\right)$。

接下来，我们对上述离散系统选取适当参数来模拟其动力学行为，设置 ${\delta }_1,{\delta }_2,{\delta }_3$ 为实参数，p为控制参数，为了便于计算取 $m=n=1$。

图2. 系统关于参数p分岔图和时间序列图，此时 ${\delta }_1=1.38,{\delta }_2=1.2,{\delta }_3=0.0035$

图3. 系统关于参数p分岔图和时间序列图，此时 ${\delta }_1=1.85,{\delta }_2=1.2,{\delta }_3=0.0035$

如图2所示，设 ${\delta }_1=1.38,{\delta }_2=1.2,{\delta }_3=0.0035$ 时，在 $p=0$ 附近系统变得不稳定，出现跳跃现象，在 $p<-0.4$ 或 $p>0.4$ 时，系统由单倍周期运动进入到二倍周期运动。在 $-0.4<p<0.4$ 时，系统做单周期运动。随着p值增大，周期1吸引子消失周期2吸引子出现。此时右图为当 $p=0.3$ 的时间序列图。

如图3所示，设 ${\delta }_1=1.85,{\delta }_2=1.2,{\delta }_3=0.0035$ 如时，图当 $p<-0.2$ 或 $p>0.2$，随着p大于临界值的时候，系统出现由周期1吸引子到周期2吸引子，并随着p增加或减少，系统进入了二倍周期运动进而到达混沌。此时右图为当 $p=0.18$ 的时间序列图。

3. 四阶Runge-Kutta算法理论过程及仿真计算

图4. 算法实现

四阶Runge-Kutta方法算法如下：

1) 输入 $a,b$ 区间等分数n，初值 $y_0$。2) 导出 $y\left(x\right)$ 在x的n个点处近似解析值y。3) 设置 $h=\frac{b-a}{n},k=0,x=a,y=y_0$。4) 计算 $K_1=f\left(x,y\right)$，$K_2=f\left(x+\frac{h}{2},y+\frac{h}{2}{K}_1\right)$，$K_3=f\left(x+\frac{h}{2},y+\frac{h}{2}{K}_2\right)$，$K_4=f\left(x+h,y+h{K}_3\right)$，置 $y+\frac{h}{6}\left(K_1+2K_2+2K_3+K_4\right)\Rightarrow y,x+h\Rightarrow x$，输出(见图4)。

对于Duffing方程的初值问题一般采用Runge-Kutta法计算，在Matlab中具体用其命令来解决初值问题，步骤如下：

① 通过实际工程模型规律，依据相关公式和定理得出相应微分方程和对应的初始条件。

② 运用相关数值变量替换，将高阶的方程变为一阶的微分方程组，并写出相应初始值。

③ 编写M文件，调整命令数据，运行得到指定区间的数值解。

对于系统(2)而言，其线性化方程为：

$\begin{array}{l}{y^{\prime }}_1=y_2\\ {y^{\prime }}_2=-m{y}_1+3n{x}_1^2{y}_1-2\mu \alpha y_2\end{array}$ (9)

令 $A=\left(\begin{array}{cc}0& 1\\ -m& -2\mu \alpha \end{array}\right),F\left(t\right)=\left(\begin{array}{cc}0& 0\\ -3n{x}_1^2& 0\end{array}\right),Y=\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)$ 方程(9)变为：

$Y^{\prime }=\left[A+F\left(t\right)\right]\cdot Y$ (10)

此时最大的Lyapunov指数为：

$\lambda =\underset{x\to \infty }{\lim }\frac{1}{t}\log \Vert \frac{Y\left(t\right)}{Y\left(0\right)}\Vert ,\Vert Y\left(t\right)\Vert =\sqrt{y_1^2+y_2^2}$ (11)

现令 $s_i=\frac{y_i}{M},i=1,2$，$M=\Vert Y\Vert =\sqrt{y_1^2+y_2^2}$ 实际计算时，步长设为 $\Delta t$，则上式化为：

$\lambda =\underset{N\to \infty }{\lim }\frac{1}{N\Delta t}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\displaystyle {\int }_{t_{i-1}}^{t_i}\frac{M^{\prime }\left(\tau \right)}{M\left(\tau \right)}\text{d}\tau }}$ (12)

然后利用四阶Runge-Kutta算法求解方程(2)、(9)即可以得出 $\lambda$ 的值。

图5. 系统的最大Lyapunov指数

例如，我们模拟参数 $\Omega =0.7,\mu f=0.1,\mu \alpha =0.15$ 时得到系统最大Lyapunov指数、相图和时间序列图，如图5、图6所示，其中 $\lambda$ 值为正，系统存在混沌现象。

图6. 当 $\Omega =0.7,\mu f=0.1,\mu \alpha =0.15$ 时，系统相图和时间序列图

4. 周期激励的振幅对Duffing系统的影响

考虑非线性Duffing方程(2)，不妨固定外激励频率 $\Omega =0.135$，且 $m,n$ 设置相应参数。对于一对给定的参数 $\mu \alpha =k_1$ 和 $\mu f=k_2$，系统存在着多种周期吸引子和混沌吸引子，因此初始条件的取值或 $k_1,k_2$ 参数的微小变化，会使得Duffing方程(2)的近似解在不同吸引子之间相互转换；故先固定阻尼值，研究周期激励振幅 $k_2$ 对Duffing系统混沌的影响。

设阻尼参数 $k_1=0.15$，系统周期激励振幅值 $k_2$ 在区间 $\left[0.01,1.5\right]$ 内分布取值，并考虑初始状态 $x_0=-0.01,y_0=0.01$ 时，分别选取参数 $k_2=0.01,k_2=0.35,k_2=0.75,k_2=1.35$，用四阶Runge-Kutta算法，Matlab程序来进行数值模拟，分别作出系统时间序列图和相图，并计算相应的最大Lyapunov指数。由时间序列图知，在阻尼参数不变的情况下，系统由不稳定的平衡点状态快速运动到左支的稳定状态，然后又回到不稳定平衡点，并向右支的稳定状态快速运动，再回到不稳定平衡点。即当 $k_1$ 取较小值时，随着 $k_2$ 值的增加，不动点由初始平衡态逐渐向左右两个吸引集运动，运动轨迹在 $\left[-3,3\right]$ 之间来回运动。

如图7(a)、图7(b)中相图所示，我们发现当 $k_2=0.01,k_2=0.35$ 时，有两个共存的混沌吸引子，且通过计算发现当 $k_2=0.01$ 时对应的最大Lyapunov指数 $\lambda =0.798791$，即 $\lambda$ 为正值，这表明该系统是混沌的。保持其它各参数不变，缓慢增大振幅值 $k_2$，研究表明它对于系统混沌状态的影响不大，如图7(c)，当 $k_2=0.75$ 或者 $k_2=1.35$ 时，仍然有两个共存的混沌吸引子，这可以通过计算最大Lyapunov指数为正得到验证，且状态大致相同，这表明在考虑各参数不变的前提下，随着系统周期激励的振幅 $k_2$ 缓慢增大，该Duffing系统混沌状态变化不明显，这也可从轨线运动范围和由平衡态进入左右两支的稳定态的时间周期得出。

图7. (a)、(b)系统在 $k_2=0.01,k_2=0.35$ 时相图和时间序列图，(c)表示系统在 $k_2=0.01,k_2=0.35$ 时的相图

5. 阻尼项对Duffing系统的影响

现研究阻尼值 $k_1$ 变化对Duffing系统的影响。设振幅参数 $k_2=0.1$，保持初始状态和其它参数不变，分别取 $k_1=0.01,k_1=0.05,k_1=0.15,k_1=0.5$，作出系统时间序列图和相图。

由图知随着阻尼参数 $k_1$ 的增加，可观察相应的混沌状态。如图8(a)~(c)所示，当 $k_1=0.01$ 系统存在混沌吸引子，做随机不规则运动，随 $k_1$ 的缓慢增大，系统在混沌状态下逐渐出现两个共存的混沌吸引子，而当 $k_1$ 逐渐增加到0.5时，共存的混沌吸引子会分离，出现单独右支的情况，如图8(d)所示。

而将 $k_1$ 继续增大，分别取 $k_1=1.2,k_1=4.5,k_1=20$，我们发现，系统由混沌吸引子逐渐向周期吸引子跃迁，随机不规则运动变为周期运动。如图9(a)、图9(b)，当 $k_1=1.2$ 时，计算最大Lyapunov指数 $\lambda =0.417286$，表明此时仍为混沌状态，而当 $k_1=4.5$ 时， $\lambda =-0.812504$ 将不再为正，吸引子在平衡点左侧会变为一个封闭轨线。如图9(c)所示， $k_1=20$ 时，系统轨线呈现出由平衡点出发作多轨运动并回到平衡点，此时 $\lambda =-0.006227$。因此，研究表明随着阻尼项的增大，系统混沌行为将发生显著变化，混沌状态将逐渐消失。

图8. (a)、(b)、(c)、(d)表示 $k_1=0.01,k_1=0.05,k_1=0.15,k_1=0.5$ 时系统相图和时间序列图

图9. (a)、(b)、(c)分别表示当 $k_1=1.2,k_1=4.5,k_1=20$ 时系统相图和时间序列图

6. 周期激励的频率对Duffing系统的影响

现研究外激励频率 $\Omega$ 对Duffing系统的影响。当系统的外激励频率远小于固有频率，即 $\Omega \ll \alpha$ 时，系统存在多时间尺度，会表现出一些复杂的快慢效应。

(c)

图11. (a)、(b)、(c)分别表示当 $\Omega =0.1355,\Omega =0.08,\Omega =0.05$ 时相图和时间序列图

7. 结论

通过分析三个参数项对这类非线性Duffing系统混沌产生的影响，我们发现，周期激励的振幅、频率或阻尼参数的改变可能导致不同的动态行为。根据四阶Runge-Kutta算法理论及相应数值模拟等研究系统产生混沌的趋势，说明这些参数的扰动会影响其状态，如振幅的改变对系统无明显影响，这可从相关实验数据得出，而随阻尼参数的增大，系统混沌行为将发生显著变化，混沌状态将逐渐消失形成周期运动。此外，对于周期激励频率的影响而言，它也可由数值模拟的相图、时间序列图及最大Lyapunov指数等相关数据表明，随着频率的增加，系统由单独左支的混沌吸引子变为左右两支共存的混沌吸引子，而将频率减小时又会由两混沌吸引子共存状态还原到单独左支情况，然后继续变为两支共存，虽然都经过了相同状态，但前后变化不完全相同。因此选择不同参数，可使系统从有序的周期解转变为混沌，也可以由混沌转变为有序态。一般来说选取较小的扰动变化来实现混沌控制，这在实际工程模型中应用较多，使操作方便简单。在非线性系统的相关研究 [18] [19] 中，通过控制其相关参数来研究系统混沌机理具有十分重要的学术意义。

基金项目

国家自然科学基金项目(11572288)和浙江省自然科学基金项目(LY20A020003)资助。

文章引用

李 亮,钱有华. 一类含周期激励Duffing系统的数值模拟
Numerical Simulation of a Class of Duffing Systems with Periodic Excitation[J]. 动力系统与控制, 2021, 10(01): 77-90. https://doi.org/10.12677/DSC.2021.101009

参考文献