DOI: 10.12677/DSC.2020.91003, PDF, HTML, XML,  被引量    国家自然科学基金支持
作者: 胡佩倩, 钱有华^*：浙江师范大学，数学与计算机科学学院，浙江 金华
关键词: 同伦分析方法；耦合van der Pol振子；周期解；Homotopy Analysis Method； Coupled van der Pol Oscillator； Periodic Solution

文章引用：胡佩倩, 钱有华. 一类环上的三自由度耦合van der Pol方程的同伦分析方法[J]. 动力系统与控制, 2020, 9(1): 22-39. https://doi.org/10.12677/DSC.2020.91003

1. 引言

随着科学技术的发展，越来越多实际问题显现出非线性的本质，故而非线性问题早已渗透于我们身边的方方面面。对于非线性问题，我们通常转变成非线性方程来解决。而对于非线性方程的求解，不像线性问题可以用线性叠加原理来求解。非线性问题求解起来相对比较复杂，而对于此问题，国内外专家学者近年来也进行了诸多研究。对于强非线性系统的定量分析方法主要有L-P摄动法 [1] 、同伦分析方法 [2] [3] [4] [5] 、双曲函数摄动法 [6] 等。

廖世俊教授于1998年首次提出同伦分析方法，并将其应用到非线性振动的研究中 [7] [8]，文献 [9] 给出非线性单自由度振动系统的两个周期公式，该公式对系统任何的振幅都是一致有效的。廖世俊等学者 [10] 应用同伦分析方法，举例并分析研究了非牛顿磁流体流动系统的解析解。钱有华等人 [11] 对同伦分析方法进行了改造，得到一类强非线性耦合van der Pol-Duffing系统的解析近似周期解。2016年，崔继峰 [12] 利用同伦分析方法并结合多尺度思想，对受迫van der Pol-Duffing振子的周期解等性质进行了研究。2016年，许天亮等人 [13] 利用同伦分析方法研究了几类简单非线性微分方程问题，并验证了同伦分析方法在求解部分非线性微分方程的优越性。2017年，郭江明 [14] 用同伦分析方法探究了带有时滞的耦合van der Pol系统双Hopf分岔和周期解。2018年，杨孟 [15] 等学者采用同伦分析方法研究了圆盘上非牛顿幂律流体的Marangoni边界层流动传热问题，利用求解化简后的非线性问题，从而获得一系列的近似解析解。

耦合van der Pol方程也是社会诸多学者研究的一个热点问题，其在动力学性质方面取得许多显著的成果。2007年，Rompala Kevin等 [16] 基于3-DOF耦合van der Pol系统，研究系统的动态行为。2009年，Jiashi Tang等 [17] 研究了对耦合van der Pol极限环的振幅进行控制的有效方法，并分析了部分弱非线性系统，得出了控制参数同振幅之间的关系，为如何有效控制振幅提出有效建议。同年，钱有华等 [18] 巧妙利用同伦分析方法探究了参数激励的非自制薄板系统的周期解。2016年，吴俊 [19] 利用KAM理论分析证明拟合周期解存在于2-DOF耦合van der Pol方程中。2018年，石兆羽等 [20] 利用混沌振子对微弱信号敏感以及对强噪声具有良好免疫力的特性，提出基于耦合van der Pol-Duffing振子系统检测微弱信号的新方法。2019年，施添添等 [21] 通过讨论特征方程根分布情况研究含时滞的大规模van der Pol-Duffing耦合振子系统的同步和异步周期振荡、概周期运动以及混沌吸引子等非线性动力学现象。

本文研究的耦合van der Pol振子环结构如图1所示 [22] [23]。而环状结构在医学、光学、电路、机电传感器和神经网络等诸多领域有着重要的应用。2007年，Fuad，Baolong Lu和Surendra Singh [24] 研究了带调制泵的双向固态环形激光器的准同步、同步和非同步混沌振荡的机制，并用最大李亚普诺夫指数进行了计算。2010年，Fernando，Ana Paula和Carla [25] 研究了由两个(单向或双向)环形细胞通过缓冲细胞耦合而成的耦合细胞网络的动力学特性。2016年，Faralli [26] 等人论证了在具有母线和环形拓扑结构的光子集成片上网络中，在同一微环中添加或删除两个波长相同的光信号的可行性。2019年，Xu，Yu和Wang [27] 为了避免简单环拓扑结构的缺点，提出了一种具有时延的耦合双环网络。采用分治算法分析相邻矩阵特征值的分布研究了周期动力学平衡和层间同步的稳定性，并进一步分析了耦合环与反射波之间的动力学相互作用。

本文主要建立三自由度耦合van der Pol振子环的精确近似解析解，且利用此系统的周期解证明了同伦分析方法的有效性和巨大的应用潜力。此外，还对同伦分析方法和数值积分法的计算结果进行了比较。结果表明，在时间历程图中，即使时间t的值在很大的区间内变化，同伦分析方法的解析解与数值积分解也能很好地吻合。

2. 同伦分析方法简述

我们先引入多自由度耦合振子的一般表达式：

$M\ddot{q}+G\stackrel{\dot{}}{q}+Kq=F\left(q\right),$ (2-1)

其中q为n维未知向量，即 $q={\left(x_1,x_2,\cdots ,x_n\right)}^{\text{T}}$，且 $\ddot{q}$ 为q关于时间变量t的二阶微分，类似地， $\stackrel{\dot{}}{q}$ 为q关于时间变量t的一阶微分，M，G，K为n阶方阵，而耦合项 $F\left(q\right)$ 则是关于变量 $\stackrel{\dot{}}{q},q,t$ 的向量值函数。

首先，我们先引入一个非线性算子：

$N\left[u\left(r,t\right)\right]=M\frac{{\partial }^{2}u\left(r,t\right)}{\partial t^2}+G\frac{\partial u\left(r,t\right)}{\partial t}+Ku\left(r,t\right)-F\left[u\left(r,t\right)\right],$ (2-2)

其中

$u\left(r,t\right)={\left(x_1\left(t\right),\cdots ,x_n\left(t\right)\right)}^{\text{T}},$ (2-3)

$\frac{\partial u\left(r,t\right)}{\partial t}={\left(\frac{\text{d}{x}_1}{\text{d}t},\cdots ,\frac{\text{d}{x}_n}{\text{d}t}\right)}^{\text{T}},$ (2-4)

$\frac{{\partial }^{2}u\left(r,t\right)}{\partial t^2}={\left(\frac{{\text{d}}^2{x}_1}{\text{d}{t}^2},\cdots ,\frac{{\text{d}}^2{x}_n}{\text{d}{t}^2}\right)}^{\text{T}},$ (2-5)

而这里的 $u\left(r,t\right)$ 是未知向量函数，r和t分别是空间与时间变量。

由同伦分析方法可知，我们需要构造一零阶形变方程

$\left(1-p\right)\left\{L\left[\Phi \left(r,t,p\right)-u_0\left(r,t\right)\right]\right\}=phH\left(t\right)N\left[\Phi \left(r,t,p\right)\right],$ (2-6)

其中 $p\in \left[0,1\right]$ 为嵌入变量，L为辅助线性算子， $u_0\left(r,t\right)$ 为初始猜测解，而h和 $H\left(t\right)$ 分别为辅助参数和辅助函数。辅助参数h对级数解的收敛性有决定性作用，故我们可画出h-ω曲线来确定h的值来确保级数解收敛。

易知，当 $p=0$ 时，零阶形变方程(2-6)演变为 $\Phi \left(r,t,0\right)=u_0\left(r,t\right)$；而当 $p=1$ 时，(2-6)又可化简为 $\Phi \left(r,t,1\right)=u\left(r,t\right)$。所以我们可得出结论：当嵌入变量p由 $0\to 1$ 的变化过程中，解 $\Phi \left(r,t,p\right)$ 从初始猜测解 $u_0\left(r,t\right)$ 变化到精确解 $u\left(r,t\right)$。

记

$u_m\left(r,t\right)=\frac{1}{m!}{\frac{{\partial }^m\Phi \left(r,t,p\right)}{\partial p^m}|}_{p=0}$ (2-7)

我们把 $\Phi \left(r,t,p\right)$ 展开成关于嵌入变量p的Taylor级数，可得

$\Phi \left(r,t,p\right)=u_0\left(r,t\right)+{\displaystyle {\sum }_{m=1}^{+\infty }{u}_m\left(r,t\right)p^m}.$ (2-8)

若选取适当的辅助线性算子L，辅助参数h和辅助函数 $H\left(t\right)$，再结合有效的初始猜测解便可使得幂级数(2-8)在 $p=1$ 时收敛，因而可得

$u\left(r,t\right)=u_0\left(r,t\right)+{\displaystyle {\sum }_{m=1}^{+\infty }{u}_m\left(r,t\right)}.$ (2-9)

向量 $u_m$ 定义为

$u_m=\left\{u_0\left(r,t\right),u_1\left(r,t\right),\cdots ,u_m\left(r,t\right)\right\}.$ (2-10)

我们对零阶形变方程(2-6)关于嵌入变量p求m次导数，再将求导得到的式子左右两边同除以 $m!$，并令 $p=0$，则我们就可得到m阶形变方程

$L\left[u_m\left(r,t\right)-{\chi }_m{u}_{m-1}\left(r,t\right)\right]=hH\left(t\right)R_m\left(u_{m-1},r,t\right),$ (2-11)

其中

${\chi }_m=\{\begin{array}{l}0,m\le 1\\ 1,m>1\end{array}$ (2-12)

且

$R_m\left(u_{m-1},r,t\right)=\frac{1}{\left(m-1\right)!}{\frac{{\partial }^{m-1}N\left(\Phi \left(r,t,p\right)\right)}{\partial p^{m-1}}|}_{p=0}.$ (2-13)

原则上m阶形变方程是线性的，而对于线性问题我们可借助数学软件如Mathematica求解。

3. 同伦分析方法的应用

本论文所研究的系统是三自由度耦合van der Pol方程，其一般表达式为

$\begin{array}{l}{\mu }_1{\ddot{x}}_1+{\epsilon }_1\left(x_1^2-\lambda \right){\stackrel{\dot{}}{x}}_1+{\eta }_1{x}_1=f_1\left(x_1,x_2,x_3\right)\\ {\mu }_2{\ddot{x}}_2+{\epsilon }_2\left(x_2^2-\lambda \right){\stackrel{\dot{}}{x}}_2+{\eta }_2{x}_2=f_2\left(x_1,x_2,x_3\right)\\ {\mu }_3{\ddot{x}}_3+{\epsilon }_3\left(x_3^2-\lambda \right){\stackrel{\dot{}}{x}}_3+{\eta }_3{x}_3=f_3\left(x_1,x_2,x_3\right)\end{array}$ (3-1)

其中 $f_1,f_2,f_3$ 互为耦合项。对于上述所给的一般方程，我们不妨规定 ${\mu }_1={\mu }_2={\mu }_3$， ${\epsilon }_1={\epsilon }_2={\epsilon }_3$， ${\eta }_1={\eta }_2={\eta }_3$，而对于方程(3-1)中 $f_i\left(i=1,2,3\right)$ 的取值，我们将根据振子运动模式分成四类来讨论 [28] ：第一，所有振子都同步运动；第二，两个振子同步运动，而第三个振子以一无关的方式运动(除它与第二个振子有相同周期的振动外)；第三，环上相邻的振子之间彼此都相位差1/3周期的运动；第四，两个振子相位差1/2周期，而第三个振子2倍于它们的频率振动。我们将针对上述的四种振子不同的运动状态研究系统的周期解问题。在文章最后，我们会采用数值模拟的形式对其进行验证，表明改进的同伦分析方法作用于此类系统是可取的，并希望该方法在其他科学领域能得到广泛应用。

要解决上述问题，我们先考虑如下的三自由度非线性耦合van der Pol振子：

$\begin{array}{l}{\mu }_1{\ddot{x}}_1+{\epsilon }_1\left(x_1^2-\lambda \right){\stackrel{\dot{}}{x}}_1+{\eta }_1{x}_1=f_1\left(x_1,x_2,x_3\right)\\ {\mu }_2{\ddot{x}}_2+{\epsilon }_2\left(x_2^2-\lambda \right){\stackrel{\dot{}}{x}}_2+{\eta }_2{x}_2=f_2\left(x_1,x_2,x_3\right)\\ {\mu }_3{\ddot{x}}_3+{\epsilon }_3\left(x_3^2-\lambda \right){\stackrel{\dot{}}{x}}_3+{\eta }_3{x}_3=f_3\left(x_1,x_2,x_3\right)\end{array}$ (3-2)

其中 ${\ddot{x}}_i\left(i=1,2,3\right)$ 为关于变量t的二阶微分， ${\stackrel{\dot{}}{x}}_i\left(i=1,2,3\right)$ 为关于变量t的一阶微分， $x_i\left(i=1,2,3\right)$ 为关于变量t的未知实函数， $f_i\left(i=1,2,3\right)$ 是耦合项 ${\mu }_i,{\epsilon }_i,{\eta }_i,\lambda \left(i=1,2,3\right)$ 为已知实参数。

我们设该系统的周期解频率为 $\omega$，然后引入一个新变量 $\tau =\omega t$，将该新变量代入系统(3-2)，得到

$\begin{array}{l}{\mu }_1{\omega }^2{\ddot{u}}_1+{\epsilon }_1\omega \left(u_1^2-\lambda \right){\stackrel{\dot{}}{u}}_1+{\eta }_1{u}_1=f_1\left(u_1,u_2,u_3\right)\\ {\mu }_2{\omega }^2{\ddot{u}}_2+{\epsilon }_2\omega \left(u_2^2-\lambda \right){\stackrel{\dot{}}{u}}_2+{\eta }_2{u}_2=f_2\left(u_1,u_2,u_3\right)\\ {\mu }_3{\omega }^2{\ddot{u}}_3+{\epsilon }_3\omega \left(u_3^2-\lambda \right){\stackrel{\dot{}}{u}}_3+{\eta }_3{u}_3=f_3\left(u_1,u_2,u_3\right)\end{array}$ (3-3)

方程满足初始条件

$u_1\left(0\right)=a,{u^{\prime }}_1\left(0\right)=b,u_2\left(0\right)=c,{u^{\prime }}_2\left(0\right)=d,u_3\left(0\right)=e,{u^{\prime }}_3\left(0\right)=0$ (3-4)

其中导数表示对变量 $\tau$ 微分。

若系统(3-3)的周期解可以用基函数 $\left\{\cos \left(k\tau \right),\sin \left(k\tau \right)|k=0,1,2,\cdots \right\}$ 表示，则有

$\begin{array}{l}{u}_1\left(\tau \right)={\displaystyle {\sum }_{k=0}^{+\infty }\left(a_{1,k}\cos k\tau +b_{1,k}\sin k\tau \right)}\\ u_2\left(\tau \right)={\displaystyle {\sum }_{k=0}^{+\infty }\left(a_{2,k}\cos k\tau +b_{2,k}\sin k\tau \right)}\\ u_3\left(\tau \right)={\displaystyle {\sum }_{k=0}^{+\infty }\left(a_{3,k}\cos k\tau +b_{3,k}\sin k\tau \right)}\end{array}$ (3-5)

我们可以假定初始猜测解可表示为以下形式

$\begin{array}{l}{u}_{1,0}\left(\tau \right)=a_0\cos \tau +b_0\sin \tau ,\\ u_{2,0}\left(\tau \right)=c_0\cos \tau +d_0\sin \tau ,\\ u_{3,0}\left(\tau \right)=e_0\cos \tau .\end{array}$ (3-6)

之后，我们还需要引入辅助线性算子

$L\left(\begin{array}{c}{\Phi }_1\left(\tau ,p\right)\\ {\Phi }_2\left(\tau ,p\right)\\ {\Phi }_3\left(\tau ,p\right)\end{array}\right)={\omega }_0^2\left(\begin{array}{c}\frac{{\partial }^2{\Phi }_1\left(\tau ,p\right)}{\partial {\tau }^2}+{\Phi }_1\left(\tau ,p\right)\\ \frac{{\partial }^2{\Phi }_2\left(\tau ,p\right)}{\partial {\tau }^2}+{\Phi }_2\left(\tau ,p\right)\\ \frac{{\partial }^2{\Phi }_3\left(\tau ,p\right)}{\partial {\tau }^2}+{\Phi }_3\left(\tau ,p\right)\end{array}\right),$ (3-7)

线性算子L满足下列性质

$L\left(\begin{array}{c}{C}_1\sin \tau +C_2\cos \tau \\ C_3\sin \tau +C_4\cos \tau \\ C_5\sin \tau +C_6\cos \tau \end{array}\right)=0.$ (3-8)

对于系统(3-2)，我们再定义非线性算子

$\begin{array}{l}N\left(\begin{array}{c}{\Phi }_1\left(\tau ,p\right)\\ {\Phi }_2\left(\tau ,p\right)\\ {\Phi }_3\left(\tau ,p\right)\end{array}\right)\\ =\left(\begin{array}{c}{\mu }_1{\Omega }^2\left(p\right)\frac{{\partial }^2{\Phi }_1\left(\tau ,p\right)}{\partial {\tau }^2}+{\epsilon }_1\Omega \left(p\right)\left[{\Phi }_1^2\left(\tau ,p\right)-\lambda \right]\frac{\partial {\Phi }_1\left(\tau ,p\right)}{\partial \tau }+{\eta }_1{\Phi }_1\left(\tau ,p\right)-f_1\left({\Phi }_1\left(\tau ,p\right),{\Phi }_2\left(\tau ,p\right),{\Phi }_3\left(\tau ,p\right)\right)\\ {\mu }_2{\Omega }^2\left(p\right)\frac{{\partial }^2{\Phi }_2\left(\tau ,p\right)}{\partial {\tau }^2}+{\epsilon }_2\Omega \left(p\right)\left[{\Phi }_2^2\left(\tau ,p\right)-\lambda \right]\frac{\partial {\Phi }_2\left(\tau ,p\right)}{\partial \tau }+{\eta }_2{\Phi }_2\left(\tau ,p\right)-f_2\left({\Phi }_1\left(\tau ,p\right),{\Phi }_2\left(\tau ,p\right),{\Phi }_3\left(\tau ,p\right)\right)\\ {\mu }_3{\Omega }^2\left(p\right)\frac{{\partial }^2{\Phi }_3\left(\tau ,p\right)}{\partial {\tau }^2}+{\epsilon }_3\Omega \left(p\right)\left[{\Phi }_3^2\left(\tau ,p\right)-\lambda \right]\frac{\partial {\Phi }_3\left(\tau ,p\right)}{\partial \tau }+{\eta }_3{\Phi }_3\left(\tau ,p\right)-f_3\left({\Phi }_1\left(\tau ,p\right),{\Phi }_2\left(\tau ,p\right),{\Phi }_3\left(\tau ,p\right)\right)\end{array}\right)\end{array}$

不妨令 $H\left(t\right)=1$，则可得零阶形变方程为

$\left(1-p\right)L\left(\begin{array}{c}{\Phi }_1\left(\tau ,p\right)-u_{1,0}\left(\tau \right)\\ {\Phi }_2\left(\tau ,p\right)-u_{2,0}\left(\tau \right)\\ {\Phi }_3\left(\tau ,p\right)-u_{3,0}\left(\tau \right)\end{array}\right)=phN\left(\begin{array}{c}{\Phi }_1\left(\tau ,p\right)\\ {\Phi }_2\left(\tau ,p\right)\\ {\Phi }_3\left(\tau ,p\right)\end{array}\right),$ (3-9)

上式满足初始条件

$\begin{array}{l}{\Phi }_1\left(0,p\right)=a\left(p\right),{\frac{\partial {\Phi }_1\left(\tau ,p\right)}{\partial \tau }|}_{\tau =0}=b\left(p\right),\\ {\Phi }_2\left(0,p\right)=c\left(p\right),{\frac{\partial {\Phi }_2\left(\tau ,p\right)}{\partial \tau }|}_{\tau =0}=d\left(p\right),\\ {\Phi }_3\left(0,p\right)=e\left(p\right),{\frac{\partial {\Phi }_3\left(\tau ,p\right)}{\partial \tau }|}_{\tau =0}=0.\end{array}$ (3-10)

显然，我们很容易得到当 $p=0$ 时，方程(3-9)满足初始条件(3-10)的解为

${\Phi }_1\left(\tau ,0\right)=u_{1,0}\left(\tau \right),{\Phi }_2\left(\tau ,0\right)=u_{2,0}\left(\tau \right),{\Phi }_3\left(\tau ,0\right)=u_{3,0}\left(\tau \right),\Omega \left(0\right)={\omega }_0.$ (3-11)

而当 $p=1$ 时，满足初始条件(3-10)的零阶形变方程(3-9)与满足初始条件(3-4)的系统(3-3)是等价的，因而有

${\Phi }_1\left(\tau ,1\right)=u_1\left(\tau \right),{\Phi }_2\left(\tau ,1\right)=u_2\left(\tau \right),{\Phi }_3\left(\tau ,1\right)=u_3\left(\tau \right),\Omega \left(1\right)=\omega .$ (3-12)

由此可见当嵌入变量p从0变化到1时， ${\Phi }_i\left(\tau ,p\right)$ 从初始猜测解 $u_{i,0}\left(\tau \right)$ 变化到未知的精确解 $u_i\left(\tau \right)\left(i=1,2,3\right)$。类似地，我们也可以得出 $\Omega \left(p\right)$ 从初始猜测频率 ${\omega }_0$ 往物理变化频率 $\omega$ 变化的结论。

我们把 ${\Phi }_1\left(\tau ,p\right),{\Phi }_2\left(\tau ,p\right),{\Phi }_3\left(\tau ,p\right)$ 与 $\Omega \left(p\right)$ 都展开成关于p的麦克劳林级数，可得到

$\begin{array}{l}{\Phi }_1\left(\tau ,p\right)=u_{1,0}\left(\tau \right)+{\displaystyle {\sum }_{m=1}^{\infty }{u}_{1,m}\left(\tau \right)p^m},\\ {\Phi }_2\left(\tau ,p\right)=u_{2,0}\left(\tau \right)+{\displaystyle {\sum }_{m=1}^{\infty }{u}_{2,m}\left(\tau \right)p^m},\\ {\Phi }_3\left(\tau ,p\right)=u_{3,0}\left(\tau \right)+{\displaystyle {\sum }_{m=1}^{\infty }{u}_{3,m}\left(\tau \right)p^m},\\ \Omega \left(p\right)={\omega }_0+{\displaystyle {\sum }_{m=1}^{\infty }{\omega }_m{p}^m}\end{array}$ (3-13)

其中

$\begin{array}{l}{u}_{1,m}\left(\tau \right)=\frac{1}{m!}{\frac{{\partial }^m{\Phi }_1\left(\tau ,p\right)}{\partial p^m}|}_{p=0},\\ u_{2,m}\left(\tau \right)=\frac{1}{m!}{\frac{{\partial }^m{\Phi }_2\left(\tau ,p\right)}{\partial p^m}|}_{p=0},\\ u_{3,m}\left(\tau \right)=\frac{1}{m!}{\frac{{\partial }^m{\Phi }_3\left(\tau ,p\right)}{\partial p^m}|}_{p=0},\\ {\omega }_m=\frac{1}{m!}{\frac{{\partial }^m\Omega \left(p\right)}{\partial p^m}|}_{p=0}.\end{array}$ (3-14)

通过选取适当的h，我们可以使(3-13)式中的幂级数在 $p=1$ 时收敛。那么，结合(3-12)式与(3-13)式，我们可以得到方程的级数解为

$\begin{array}{l}{u}_1\left(\tau \right)=u_{1,0}\left(\tau \right)+{\displaystyle {\sum }_{m=1}^{\infty }{u}_{1,m}}\left(\tau \right),\\ u_2\left(\tau \right)=u_{2,0}\left(\tau \right)+{\displaystyle {\sum }_{m=1}^{\infty }{u}_{2,m}}\left(\tau \right),\\ u_3\left(\tau \right)=u_{3,0}\left(\tau \right)+{\displaystyle {\sum }_{m=1}^{\infty }{u}_{3,m}}\left(\tau \right),\\ \omega ={\omega }_0+{\displaystyle {\sum }_{m=1}^{\infty }{\omega }_m}.\end{array}$ (3-15)

记

$\begin{array}{l}{u}_{1,n}=\left\{u_{1,0}\left(\tau \right),u_{1,1}\left(\tau \right),\cdots ,u_{1,n}\left(\tau \right)\right\},\\ u_{2,n}=\left\{u_{2,0}\left(\tau \right),u_{2,1}\left(\tau \right),\cdots ,u_{2,n}\left(\tau \right)\right\},\\ u_{3,n}=\left\{u_{3,0}\left(\tau \right),u_{3,1}\left(\tau \right),\cdots ,u_{3,n}\left(\tau \right)\right\},\\ {\omega }_n=\left\{{\omega }_0,{\omega }_1,\cdots ,{\omega }_n\right\}.\end{array}$ (3-16)

同样地，我们对零阶形变方程(3-9)关于嵌入变量p求m次导数，再将求导得到的式子左右两边同除以 $m!$，并令 $p=0$，则我们就可得到m阶形变方程

$L\left(\begin{array}{c}{u}_{1,m}\left(\tau \right)-{\chi }_m{u}_{1,m-1}\left(\tau \right)\\ u_{2,m}\left(\tau \right)-{\chi }_m{u}_{2,m-1}\left(\tau \right)\\ u_{3,m}\left(\tau \right)-{\chi }_m{u}_{3,m-1}\left(\tau \right)\end{array}\right)=h\left(\begin{array}{c}{R}_{1,m}\left(u_{1,m-1},{\omega }_{m-1}\right)\\ R_{2,m}\left(u_{2,m-1},{\omega }_{m-1}\right)\\ R_{3,m}\left(u_{3,m-1},{\omega }_{m-1}\right)\end{array}\right),$ (3-17)

其中

$\left(\begin{array}{c}{R}_{1,m}\left(u_{1,m-1},{\omega }_{m-1}\right)\\ R_{2,m}\left(u_{2,m-1},{\omega }_{m-1}\right)\\ R_{3,m}\left(u_{3,m-1},{\omega }_{m-1}\right)\end{array}\right)=\frac{1}{\left(m-1\right)!}{\frac{{\partial }^{m-1}N\left(\begin{array}{c}{\Phi }_1\left(\tau ,p\right)\\ {\Phi }_2\left(\tau ,p\right)\\ {\Phi }_3\left(\tau ,p\right)\end{array}\right)}{\partial p^{m-1}}|}_{p=0},$ (3-18)

此m阶形变方程满足以下初始条件

$\begin{array}{l}{u}_{1,m}\left(0\right)=a_m,{u^{\prime }}_{1,m}\left(0\right)=b_m,u_{2,m}\left(0\right)=c_m,\\ {u^{\prime }}_{2,m}\left(0\right)=d_m,u_{3,m}\left(0\right)=e_m,{u^{\prime }}_{3,m}\left(0\right)=0,\left(m\ge 1\right)\end{array}$ (3-19)

由线性算子L的定义与解的表达原则可知，方程(3-17)的等式右边项里不能含有 $\cos \tau$ 与 $\sin \tau$，也不能含有久期项 $\tau \cos \tau$ 与 $\tau \sin \tau$。

我们把这些项都找出来，且令这些项的系数都为0，则可得

$\begin{array}{l}\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{1,m}\left(u_{1,m-1},{\omega }_{m-1}\right)\cos \tau \text{d}\tau }=0,\\ \frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{1,m}\left(u_{1,m-1},{\omega }_{m-1}\right)\sin \tau \text{d}\tau }=0,\\ \frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{2,m}\left(u_{2,m-1},{\omega }_{m-1}\right)\cos \tau \text{d}\tau }=0,\\ \frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{2,m}\left(u_{2,m-1},{\omega }_{m-1}\right)\sin \tau \text{d}\tau }=0,\\ \frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{3,m}\left(u_{3,m-1},{\omega }_{m-1}\right)\cos \tau \text{d}\tau }=0,\\ \frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{3,m}\left(u_{3,m-1},{\omega }_{m-1}\right)\sin \tau \text{d}\tau }=0,\end{array}$ (3-20)

其中

$\begin{array}{l}{R}_{1,m}\left(u_{1,m-1},{\omega }_{m-1}\right)={\displaystyle {\sum }_{n=1}^{{\phi }_1\left(m\right)}\left[a_{m,n}^1\cos n\tau +b_{m,n}^1\sin n\tau \right]},\\ R_{2,m}\left(u_{2,m-1},{\omega }_{m-1}\right)={\displaystyle {\sum }_{n=1}^{{\phi }_2\left(m\right)}\left[a_{m,n}^2\cos n\tau +b_{m,n}^2\sin n\tau \right]},\\ R_{3,m}\left(u_{3,m-1},{\omega }_{m-1}\right)={\displaystyle {\sum }_{n=1}^{{\phi }_2\left(m\right)}\left[a_{m,n}^3\cos n\tau +b_{m,n}^3\sin n\tau \right]}.\end{array}$ (3-21)

且

$\begin{array}{l}{a}_{m,n}^1=\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{1,m}\left(u_{1,m-1},{\omega }_{m-1}\right)\sin \left(n\tau \right)\text{d}\tau },\\ b_{m,n}^1=\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{1,m}\left(u_{1,m-1},{\omega }_{m-1}\right)\sin \left(n\tau \right)\text{d}\tau },\\ a_{m,n}^2=\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{2,m}\left(u_{2,m-1},{\omega }_{m-1}\right)\cos \left(n\tau \right)\text{d}\tau },\\ b_{m,n}^2=\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{2,m}\left(u_{2,m-1},{\omega }_{m-1}\right)\sin \left(n\tau \right)\text{d}\tau },\\ a_{m,n}^3=\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{3,m}\left(u_{3,m-1},{\omega }_{m-1}\right)\cos \left(n\tau \right)\text{d}\tau },\\ b_{m,n}^3=\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{R}_{3,m}\left(u_{3,m-1},{\omega }_{m-1}\right)\sin \left(n\tau \right)\text{d}\tau },\end{array}$ (3-22)

根据(3-17)和(3-20)式，我们可以求出 ${\omega }_k,a_k,b_k,c_k,d_k,e_k\left(k=0,1,2,\cdots \right)$ 的值。为使解的精度更高，我们需要对w进行一定程度的修正

${\omega }_{-}=\omega +v.$

其中v为小参数。

4. 数值模拟与讨论

为了验证同伦分析方法应用于非线性多自由度系统的正确性和有效性，针对耦合van der Pol系统(3-2)，我们将分四种情形对其数值模拟进行讨论，并对分别对这四种情形下的理论结果与数值结果进行比较。而耦合项 $f_i\left(i=1,2,3\right)$ 的取值可借鉴文献 [29] [30]。

4.1. 所有振子都同步运动

令

$\begin{array}{l}{f}_1\left(x_1,x_2,x_3\right)=x_1-x_2\\ f_2\left(x_1,x_2,x_3\right)=x_2-x_3\\ f_3\left(x_1,x_2,x_3\right)=x_3-x_1\end{array}$ (4-1)

我们取 ${\mu }_i=1,{\epsilon }_i=1,{\eta }_i=1,\lambda =1\left(i=1,2,3\right)$，且令 $a,b,c,d,e,\omega$ 的初始猜测值分别为 $a_0=2,b_0=0,c_0=2,d_0=2,e_0=2,{\omega }_0=1$，我们需要选取合适的收敛控制参数h来确保级数解的收敛。为此，我们给出了图2所示的ω-h曲线来确定h的取值范围。细观图2可知，当 $-1.85<h<-0.25$ 时，级数解 $\omega$ 是收敛的。为简便起见，我们取 $h=-1$，图3和图4分别给出了八阶近似相图和时间历程曲线。

利用相图(图3)与时间历程图(图4)可证明同伦分析方法求得的相图与数值分析方法所求得的相图结果比较吻合，且能明显看出 $x_1\left(\tau \right),x_2\left(\tau \right),x_3\left(\tau \right)$ 三个振子同步运动。其中初始条件为：

$\begin{array}{l}{x}_1\left(0\right)=2.008600,{x^{\prime }}_1\left(0\right)=0,\\ x_2\left(0\right)=2.008600,{x^{\prime }}_2\left(0\right)=0,\\ x_3\left(0\right)=2.008600,{x^{\prime }}_3\left(0\right)=0,\end{array}$

根据文献 [31] 并且可求得 $x_1\left(\tau \right),x_2\left(\tau \right),x_3\left(\tau \right),\omega$ 的八阶近似表达式为

$\begin{array}{c}{x}_1\left(\tau \right)=1.881842\cos \bar{\omega }\tau +0.720593\sin \bar{\omega }\tau +0.720593\cos 3\bar{\omega }\tau -1.881842\sin 3\bar{\omega }\tau \\ +0.158456\cos 5\bar{\omega }\tau -0.530028\sin 5\bar{\omega }\tau -0.176676\cos 7\bar{\omega }\tau -0.475368\sin 7\bar{\omega }\tau \\ -0.020655\cos 9\bar{\omega }\tau -0.218080\sin 9\bar{\omega }\tau -0.043616\cos 11\bar{\omega }\tau +0.103273\sin 11\bar{\omega }\tau \\ -0.010428\cos 13\bar{\omega }\tau -0.000519\sin 13\bar{\omega }\tau -0.000074\cos 15\bar{\omega }\tau +0.072994\sin 15\bar{\omega }\tau ,\end{array}$

$\begin{array}{c}{x}_2\left(\tau \right)=1.881842\cos \bar{\omega }\tau +0.720593\sin \bar{\omega }\tau +0.720593\cos 3\bar{\omega }\tau -1.881842\sin 3\bar{\omega }\tau \\ +0.158456\cos 5\bar{\omega }\tau -0.530028\sin 5\bar{\omega }\tau -0.176676\cos 7\bar{\omega }\tau -0.475368\sin 7\bar{\omega }\tau \\ -0.020655\cos 9\bar{\omega }\tau -0.218080\sin 9\bar{\omega }\tau -0.043616\cos 11\bar{\omega }\tau +0.103273\sin 11\bar{\omega }\tau \\ -0.010428\cos 13\bar{\omega }\tau -0.000519\sin 13\bar{\omega }\tau -0.000074\cos 15\bar{\omega }\tau +0.072994\sin 15\bar{\omega }\tau ,\end{array}$

$\begin{array}{c}{x}_3\left(\tau \right)=1.881842\cos \bar{\omega }\tau +0.720593\sin \bar{\omega }\tau +0.720593\cos 3\bar{\omega }\tau -1.881842\sin 3\bar{\omega }\tau \\ +0.158456\cos 5\bar{\omega }\tau -0.530028\sin 5\bar{\omega }\tau -0.176676\cos 7\bar{\omega }\tau -0.475368\sin 7\bar{\omega }\tau \\ -0.020655\cos 9\bar{\omega }\tau -0.218080\sin 9\bar{\omega }\tau -0.043616\cos 11\bar{\omega }\tau +0.103273\sin 11\bar{\omega }\tau \\ -0.010428\cos 13\bar{\omega }\tau -0.000519\sin 13\bar{\omega }\tau -0.000074\cos 15\bar{\omega }\tau +0.072994\sin 15\bar{\omega }\tau ,\end{array}$

$\omega =\mathrm{0.942953.}$

4.2. 两个振子同步运动，而第三个振子以一无关的方式运动(除它与第二个振子有相同周期地振动外)

令

$\begin{array}{l}{f}_1\left(x_1,x_2,x_3\right)=x_1-x_2\\ f_2\left(x_1,x_2,x_3\right)=x_2-x_1\\ f_3\left(x_1,x_2,x_3\right)=x_1-x_2-10x_3\end{array}$ (4-2)

在此情况中，我们取和情形4.1一样的参数值，即 ${\mu }_i=1,{\epsilon }_i=1,{\eta }_i=1,\lambda =1\left(i=1,2,3\right)$， $a,b,c,d,e,\omega$ 的初始猜测值取为 $a_0=2,b_0=0,c_0=2,d_0=0,e_0=2,{\omega }_0=1$，我们可以从图5中确定出h的取值范围来确保级数解的收敛。因为h的取值范围为 $-1.05<h<-0.65$，为简便起见故仍然取 $h=-1$，一阶近似相图与时间历程图如图6与图7所示。

利用相图(图6)与时间历程图(图7)可证明同伦分析方法求得的相图与数值分析方法所求得的相图结果比较吻合。观察图7可发现，振子 $x_1\left(\tau \right)$ 与 $x_2\left(\tau \right)$ 同步运动，而第三个振子 $x_3\left(\tau \right)$ 则以一种无关于前两个振子的方式运动。其初始条件为 $x_1\left(0\right)=4,{x^{\prime }}_1\left(0\right)=0,x_2\left(0\right)=4,{x^{\prime }}_2\left(0\right)=0,x_3\left(0\right)=4,{x^{\prime }}_3\left(0\right)=0$，且可得 $x_1\left(\tau \right),x_2\left(\tau \right),x_3\left(\tau \right),\omega$ 近似表达式为：

$\begin{array}{c}{x}_1\left(\tau \right)=1.844822\cos \bar{\omega }\tau -0.118627\sin \bar{\omega }\tau -0.273745\cos 3\bar{\omega }\tau -1.844822\sin 3\bar{\omega }\tau \\ +0.453207\cos 5\bar{\omega }\tau -0.120863\sin 5\bar{\omega }\tau -1.836718\cos 7\bar{\omega }\tau -0.135962\sin 7\bar{\omega }\tau \\ -0.047096\cos 9\bar{\omega }\tau -0.125599\sin 9\bar{\omega }\tau -0.144917\cos 11\bar{\omega }\tau +0.235481\sin 11\bar{\omega }\tau \\ +0.500366\cos 13\bar{\omega }\tau -0.003442\sin 13\bar{\omega }\tau -0.021553\cos 15\bar{\omega }\tau -0.035026\sin 15\bar{\omega }\tau ,\end{array}$

$\begin{array}{c}{x}_2\left(\tau \right)=1.844822\cos \bar{\omega }\tau -0.118627\sin \bar{\omega }\tau -0.273745\cos 3\bar{\omega }\tau -1.844822\sin 3\bar{\omega }\tau \\ +0.453207\cos 5\bar{\omega }\tau -0.120863\sin 5\bar{\omega }\tau -1.836718\cos 7\bar{\omega }\tau -0.135962\sin 7\bar{\omega }\tau \\ -0.047096\cos 9\bar{\omega }\tau -0.125599\sin 9\bar{\omega }\tau -0.144917\cos 11\bar{\omega }\tau +0.235481\sin 11\bar{\omega }\tau \\ +0.500366\cos 13\bar{\omega }\tau -0.003442\sin 13\bar{\omega }\tau -0.021553\cos 15\bar{\omega }\tau -0.035026\sin 15\bar{\omega }\tau ,\end{array}$

$\begin{array}{c}{x}_3\left(\tau \right)=1.124328\cos \bar{\omega }\tau +0.573061\sin \bar{\omega }\tau +0.134443\cos 3\bar{\omega }\tau -1.131251\sin 3\bar{\omega }\tau \\ +0.359854\cos 5\bar{\omega }\tau -0.670447\sin 5\bar{\omega }\tau -0.001711\cos 7\bar{\omega }\tau -0.223488\sin 7\bar{\omega }\tau \\ +0.271384\cos 9\bar{\omega }\tau +0.348288\sin 9\bar{\omega }\tau +0.054345\cos 11\bar{\omega }\tau -0.059881\sin 11\bar{\omega }\tau \\ -0.171917\cos 13\bar{\omega }\tau -0.066372\sin 13\bar{\omega }\tau +0.167520\cos 15\bar{\omega }\tau -0.089623\sin 15\bar{\omega }\tau ,\end{array}$

$\omega =\mathrm{1.208175.}$

4.3. 环上相邻的振子之间彼此都相位差1/3周期的运动

令

$\begin{array}{l}{f}_1\left(x_1,x_2,x_3\right)=\sqrt{3}/9{\left(x_1+2x_2\right)}^3\\ f_2\left(x_1,x_2,x_3\right)=\sqrt{3}/9{\left(x_2+2x_3\right)}^3\\ f_3\left(x_1,x_2,x_3\right)=\sqrt{3}/9{\left(x_3+2x_1\right)}^3\end{array}$ (4-3)

在此情况中，环上相邻的振子之间彼此都相位差 周期的运动，因而对系统(3-1)中的各系数取值与前两种情况有所不同。这里我们取 ${\mu }_i=1,{\epsilon }_i=1,{\eta }_i=4,\lambda =1\left(i=1,2,3\right)$。同样地，初始猜测值也被赋予一些新的数值， $a_0=-\sqrt{10}/5,b_0=\sqrt{30}/5,c_0=-\sqrt{10}/5,d_0=-\sqrt{30}/5,e_0=2\sqrt{10}/5,{\omega }_0=2$。为确保级数解的收敛，我们取 $h=-1$，从而得到相应地相图与时间历程图(见图8与图9)。

由图8可知，利用同伦分析方法求得的近似解与数值积分法求得的解基本吻合，而由图10也可看出三个振子正相差1/3周期运动。我们同样可求得数值积分法的初始条件为

$x_1\left(0\right)=-0.632455,{x^{\prime }}_1\left(0\right)=1.095445,x_2\left(0\right)=-0.632455,{x^{\prime }}_2\left(0\right)=-1.095445,x_3\left(0\right)=1.264911,{x^{\prime }}_3\left(0\right)=0$。

我们也就能得到 $x_1\left(\tau \right),x_2\left(\tau \right),x_3\left(\tau \right)$ 与 $\omega$ 的近似解表达式：

$x_1\left(\tau \right)=-0.649572\cos \tau +1.089516\sin \tau +0.017116\cos 3\tau +0.001976\sin 3\tau ,$

$x_2\left(\tau \right)=-0.615339\cos \tau -1.101374\sin \tau -0.017116\cos 3\tau +0.001976\sin 3\tau ,$

$x_3\left(\tau \right)=1.264911\cos \tau +0.094868\sin \tau -0.031623\sin 3\tau ,$

$\omega =2.$

4.4. 两个振子相位差1/2周期，而第三个振子2倍于它们的频率振动

令

$\begin{array}{l}{f}_1\left(x_1,x_2,x_3\right)=\left(1-x_1^2\right){x^{\prime }}_2\\ f_2\left(x_1,x_2,x_3\right)=\left(1-x_2^2\right){x^{\prime }}_1\\ f_3\left(x_1,x_2,x_3\right)=x_1^3{x^{\prime }}_2^3\end{array}$ (4-4)

此时，两个振子相位差1/2周期，而第三个振子 倍于它们的频率振动。为计算方便，我们取 ${\mu }_i=1,{\epsilon }_i=1,{\eta }_i=1,\lambda =1\left(i=1,2\right),{\mu }_3=1,{\epsilon }_3=1/16,{\eta }_3=4$，

而初始猜测值求得为 $a_0=1,b_0=0,c_0=-1,d_0=0,e_0=1,{\omega }_0=1$。取 $h=-1$ 得近似解的相图和时间历程图如图11与图12所示。

通过观察与计算可知，利用同伦法得到的相图与时间历程曲线是合适的，再次证明了同伦分析方法于耦合多自由度振子环系统上的适用性。其数值积分法满足的初始条件为

$x_1\left(0\right)=1,{x^{\prime }}_1\left(0\right)=0,x_2\left(0\right)=-1,{x^{\prime }}_2\left(0\right)=0,x_3\left(0\right)=1,{x^{\prime }}_3\left(0\right)=0$。

同样地，得到 $x_1\left(\tau \right),x_2\left(\tau \right),x_3\left(\tau \right)$ 和频率 $\omega$ 的一阶近似解为：

$x_1\left(\tau \right)=\cos \tau ,$

$x_2\left(\tau \right)=-\cos \tau ,$

$x_3\left(\tau \right)=\cos 2\tau ,$

$\omega =1.$

显然，此近似解就是满足条件(4-4)的系统(3-1)的精确解，也符合第四种情况的振子运动(见图13)。

5. 结论

本文利用该方法对多自由度耦合van der Pol振子环进行了精确的近似分析。由之前的数值模拟可以知道，利用同伦分析方法求解三自由度耦合van der Pol方程得到的周期解比较精确。该方法不仅可以计算任意阶次的频率和位移，而且计算量小。与经典摄动技术相比，该方法不需要依靠小参数并可以很好地保证级数解的收敛性以及解的精度，在解决非线性问题方面具有非常好的效果。在后续研究中，我们可以考虑四种不同情况中的振子如何才能互相转换，在何种外激励下使系统形成四稳态局面。

致谢

感谢国家自然科学基金项目(11572288)、浙江省自然科学基金项目(LY20A020003)和浙江师范大学数计学院重点项目(ZD01201713)资助。

NOTES

^*通讯作者。