Duffing方程的辛精细积分方法研究*

都琳侯平兰

（1.西北工业大学理学院，西安 710072）（2.陕西省建筑职工大学，西安 710068）

Duffing方程的辛精细积分方法研究*

都琳1†侯平兰2

（1.西北工业大学理学院，西安 710072）（2.陕西省建筑职工大学，西安 710068）

辛精细积分方法汲取了辛几何算法保持动力学系统辛结构的优点和精细积分方法高精度的数值优点，其实现过程中涉及到大量矩阵求逆运算.为减小辛精细积分方法的运算量，本文在辛精细积分算法之前先将非齐次方程近似齐次化，使得矩阵求逆部分不显含时间，降低矩阵求逆计算量，并将这一方法应用于无阻尼Duffing方程的数值分析.通过与经典四阶Runge-Kutta格式及精细积分方法对比，发现辛精细积分方法在数值精度、计算耗时、保持系统能量等方面明显优于Runge-Kutta格式.此外，与精细积分方法相比，辛精细积分方法在保持系统能量方面存在明显优势.

辛精细积分方法，Duffing方程，齐次化，辛几何

引言

30多年前，冯康教授等针对有限维Hamilton系统的长时间数值积分，基于辛几何理论，创立了辛算法［1］.辛算法具有长时间数值稳定性的奥妙在于其在算法设计过程中，每一步都能精确保持系统辛结构.正是由于这一开创性的研究成果，冯康教授及其研究团队荣获1997年自然科学一等奖.此后，辛算法理论框架逐渐完善［2］，并在天体力学等领域得到了一定的推广应用［3］.

在长期的研究过程中，力学工作者发现，辛算法在以下方面存在不足：一是辛几何理论太过于数学化，何谓“辛”一直是将辛算法理论应用于力学问题的瓶颈；二是辛算法本身与系统的对称性密切相关［4］，也就是说在构造辛算法之前，对系统的对称性需要有深刻的认识.基于以上两点，那么有没有一种能够适用于力学系统的时程积分方法呢？针对这一问题，我国钟万勰教授针对描述动力学问题的微分方程模型，创立了具有极高数值精度的精细积分方法［5］，这一开创性的研究工作被公认为是近年来计算力学领域最重要的研究成果之一，在动力学与控制领域、弹性力学领域得到了较为广泛的应用［6，7］，钟万勰教授研究团队也因此荣获2010年自然科学二等奖.

钟万勰教授创立的精细积分方法，回避了辛几何理论中复杂的数学理论，实用性更强了.然而，也正因为如此，精细积分方法能否保辛一直是学术界讨论的焦点，从数值结果上看，钟万勰教授认为精细积分方法是近似保辛的.为了让精细积分方法从是否保辛这一议题中解脱出来，邓子辰教授和黄永安教授等将辛算法思想引入精细积分方法，提出了辛精细积分方法［8］，使得精细积分方法在提高数值计算精度的同时，能够长时间保持系统辛结构，提高数值算法的稳定性.

对于非齐次微分方程系统，辛精细积分方法处理起来计算量将大幅增加，针对这一问题，本文首先将凑微分思想应用于非齐次微分方程系统，近似地将非齐次方程系统转化为齐次系统，然后再采用辛精细积分方法求解，得到了具有长时间数值稳定和高精度的数值分析结果，同时计算量较直接使用辛精细积分方法大幅降低.

1 受迫Duffing方程及其齐次化

钟万勰教授从一根弹簧的振动入手，深入浅出地阐述了力、功、能量之间的本质联系.本部分将在回顾力、功及能量之间本质联系的基础上，回顾Duffing方程的力学背景，并采用凑微分思想，将含有外激励的Duffing方程这一非齐次方程近似齐次化，为后续开展Duffing方程的辛精细积分方法研究奠定基础.

众所周知，弹性系统的胡克定律中，弹簧的弹性势能U与位移（形变）x的二次幂成正比，以一维情形为例，即：

其中x表示弹簧振子离开平衡位置的位移，ω为常数.式（1）的微分形式表明系统的恢复力f与位移成正比（线性关系）：

这一线性运动方程成立的前提条件是弹簧系统中的应力与应变之间成线性关系.

然而，实际的许多弹簧系统（如建筑工程中的大量构件和桁架）并不服从这样简单的规律.一般来讲，一维弹性系统的势能具有如下的普遍形式：

其中，κ，λ，μ均为常数，由此得到系统的恢复力形式为：

上式是典型的非线性恢复力.虽然实际的弹簧系统大多是非线性的，但是其弹性势能却是具有某种对称性，其中最简单的形式是只取弹性势能前两项：

通常，与κ相比，μ是一个小量，此时的恢复力为：

这种具有对称性（f的奇偶性）的弹性势能对应的运动微分方程为：

这就是著名的无阻尼Duffing方程.在周期性的外力作用下，Duffing方程将写成：

其中F，Ω分别为外力的振幅和频率.对于式（9）这类非齐次微分方程系统，辛精细积分方法的计算量较大，为此，本文将首先采用凑微分思想［9］，将式（9）近似齐次化.

在周期外力作用下的无阻尼Duffing方程（9）利用凑微分思想［9］可以变为如下的齐次方程的形式：

其中y＝x－（k1cosΩt＋k2sinΩt），k1，k2均为待定系数.将其代入式（10），得：

部分展开后，忽略部分高阶项得：

由于μ远小于κ，所以可做如下的近似（实际上，可以直接将μ［x－（k1cosΩt＋k2sinΩt）］3这一项近似为零，但是为了使用凑微分思想，保留了μx3这一可能比μ（k1cosΩt＋k2sinΩt）3这一项值还小的非线性项）：

将公式（13）代入公式（12）得到齐次方程（10）的近似形式：

与式（9）比较系数后得到待定系数的值为：

从而得到周期外力作用下的无阻尼达芬方程（9）近似等价于齐次微分方程：

即引入的新的变量可以定义为：

2 辛精细积分格式

针对方程（16）的辛精细积分方法基本思路［8］：将方程（16）降阶后，其系数矩阵可以写成如下形式：

其中，J2n为反对称阵，Hz为方阵（含变量z），H0为定常矩阵，H1（z）为非定常部分（包含非线性项和非齐次项，在这里只是包含非线性项，因为方程（16）是齐次的），z为状态向量，对于方程（16），则方程（16）可以写成：

在时间区段［tk，tk＋1］内，将非定常部分线性化处理，具体处理方法及过程参见文献［10］：

首先求解齐次线性方程：

Hamilton矢量场有一个重要性质，即对于足够小的时间间隔t，exp（H0·t）一定是一个正则映射，也就是说，exp（H0·t）应该能够精确保持哈密尔顿系统的辛结构，即exp（H0·t）为一个辛映射.

回顾精细积分法的全过程，不难发现：精细积分法的高数值精度来源于一系列的乘法变换.然而，在辛空间内，两个辛矩阵的求和计算并不能保辛，即加法不能保辛，但是两个辛矩阵的乘积结构还是辛矩阵，即乘法是保辛的.辛矩阵的这一性质为我们构造辛精细积分算法提供了一条思路，将指数矩阵函数exp（H·Δt）进行合理的辛近似，得到的指数矩阵函数exp（H·τ）就会具有辛算法的性质，能够实现系统的近似保辛.实现这一思想的关键问题在于用数值方法精确地计算变换矩阵exp（H·τ）.

利用有理Padé逼近，exp（H·τ）的Padé逼近具体形式如下［8］：

式中：

已有研究结果表明［10］：设H是无穷小的辛阵，对于充分小的，任意方阵glm（H·τ）是辛阵，当且仅当l＝m，即Ωll（x）是对Padé对角逼近的.

依据上述结论，可以通过选择l＝m＝2，3，…，就可以使得公式（21）近似保辛.同时，不难证明l，m取值越大，逼近效果越好.为了方便说明构造辛精细积分算法，选取l＝m＝2，则：

从而得到传递函数形式：

其中，Ω（H0·τ）是一个具有固定精度的数值格式，本文后续数值实验将采用四阶精度格式为例［8］：

利用精细积分思想，可以得到矩阵函数exp（H0· τ）的逼近公式：

其中

式（25）就可以方便地采用精细积分方法进行数值积分了.联合非定常项的线性近似（22），就可以得到非线性方程（16）解的迭代格式为：

需要说明的是，上述格式中，矩阵求逆部分不显含时间，这加速了矩阵求逆这一计算量较大的运算过程.

3 无阻尼Duffing方程的数值实验结果

选取参数：κ＝1，μ＝1.0×10－6，F＝0.2，Ω＝1.5，取时间步长τ＝0.2，分别用辛精细积分算法、精细积分算法（参考文献［10］中非线性非齐次方程的精细积分格式）和传统的四阶Runge-Kutta格式求解周期外力作用下无阻尼Duffing方程（9），得到弹簧系统的位移x，速度dx／dt和加速度d2x／dt2的对比情况如图1～3.

三种积分方法的计算时间见表1.从表中不难看出，精细积分方法和辛精细积分方法的计算时间远小于四阶Runge-Kutta格式的求解时间，辛精细积分方法的计算时间略小于精细积分方法的计算时间.其原因是辛精细积分方法执行之前，我们先对非齐次方程近似齐次化，降低了辛精细积分过程中矩阵求逆的计算量.

由于考虑的是无阻尼的情况，因此原系统是一个保守系统，不存在能量耗散.为了反映辛精细积分算法能够保持系统的不存在能量耗散这一特性，求解中特别记录了三种算法每一步的系统能量损失（即弹性势能与单位质量弹簧振子振动动能之和的增量）情况，如图4.

图1 三种积分方法得到的系统位移Fig.1 Displacements obtained by the three integral methods

图2 三种积分方法得到的系统速度Fig.2 Speeds obtained by the three integral methods

图3 三种积分方法得到的系统加速度Fig.3 Accelerations obtained by the three integral methods

表1 三种积分格式的计算耗时对比Table 1 Time-consuming of three integral methods

图4 三种积分方法得到的系统能量损失情况对比Fig.4 Energy loss obtained by the three integral methods

由分别利用辛精细积分算法、精细积分算法和四阶Runge-Kutta格式求解方程（9）得到系统位移、速度和加速度波形对比中不难发现：辛精细积分算法与精细积分算法计算精度相当，并且明显高于四阶Runge－Kutta格式.另外，从能量损失情况看，精细积分算法虽然具有较高的数值精度，但是它与传统的四阶Runge－Kutta格式类似，不能保持系统的几何特性，而辛精细积分算法能够较好地长时间保持系统的几何特性（能量）.

3 结论

辛算法由于其良好的长时间数值稳定性和保辛特点已被应用于力学和天文学等领域，而精细积分方法注重于算法精度和计算速度，是一种更为通俗易懂的动力学问题计算方法.结合二者的优点，辛精细积分方法将具有重要的应用价值.

本文针对辛精细积分方法实现过程中面临的由于矩阵求逆增加的计算量问题，借鉴凑微分思想将非齐次微分方程近似齐次化，使得辛精细积分格式中的非定常项系数矩阵求逆不显含时间变量，降低矩阵求逆计算量，并将这一方法应用于无阻尼Duffing方程的数值分析.数值分析结果表明：经过近似齐次化后的辛精细积分方法在数值精度、计算耗时和保持系统能量等方面明显优于经典Runge－Kutta方法；同时，在保持系统能量方面，对现有的精细积分方法有明显改进.

1 Feng K.On difference schemes and symplectic geometry，In：Proceedings of the 5th International.Symposium on Differential Geometry＆Differential Equations.Beijing：Science Press，1984：42～58

2 冯康，秦孟兆.Hamilton动力学体系的Hamilton算法.自然科学进展-国家重点实验室通讯，1991，2：102～112（Feng K，Qin M Z.Hamilton algorithm of Hamilton mechanical systems.Progress in Natural Science-State Key Laboratory of Communications，1991，2：102～112（in Chinese））

3 刘林.天体力学方法.南京：南京大学出版社，1998（Liu L.Celestial mechanics-methods.Nanjing：Nanjing University，1998（in Chinese））

4 Marsden J E，Ratiu T S.Introduction to mechanics and symmetry.New York：Springer-Verlag，1994

5 Zhong W X.On Precise integration method.Journal of Computational and Applied Mathematics，2004，163（1）：59～78

6 Zhong W X，Lin JH，Gao Q.The precise computation for wave propagation in stratified materials.International Journal for Numerical Methods in Engineering，2004，60（1）：11～25

7 Gao Q，Lin JH，Zhong W X，Howson W P，Williams F W.A precise numerical method for rayleigh waves in a stratified half space.International Journal for Numerical Methods in Engineering，2006，67（6）：771～786

8 Huang Y A，Deng ZC，Yao L X.An improved symplectic precise integration method for analysis of the rotating rigidflexible coupled system.Journal of Sound and Vibration，2007，299（1-2）：229～24

9 胡伟鹏.一类二阶常系数非齐次微分方程的解法.高等数学研究，2001，4（2）：45～46（Hu W P.Mehtod of a class of second order inhomogeneous differential equation with constant coefficients.Studies in College Mathematics，2001，4（2）：45～46（in Chinese））

10 钟万勰.应用力学的辛数学方法.北京：高等教育出版社，2006（Zhong W X.Symplectic solution methodology in applied mechanics.Beijing：Higher Education Press（in Chinese） ）

Received 27 July 2016，revised 21 August2016.

*The project supported by the National Natural Science Foundation of China（11672233，11302169）

†Corresponding author E-mail：lindu＠nwpu.edu.cn

SYMPLECTIC PRECISE INTEGRATION METHOD FOR DUFFING EQUATION*

Du Lin1†Hou Pinglan2
（1.School of Science，Northwestern Polytechnical University，Xi′an 710072，China）（2.Architecture Zabor University of Shaanxi Province，Xi′an 710068，China）

The symplectic precise integration method owns the advantages of the symplectic method and the precise integration method.In the implementation procedure of which，matrix inversion is a time-consuming step.Aiming at this problem，we homogenize the inhomogeneous equation approximately before the design of the symplectic precise integration method in this paper.The homogenizing process makes the matrix inversion time-independent and reduces the calculated quantity of the matrix inversion process，which is used in the symplectic precise integration method of the non-damping Duffing equation in this paper.From the numerical results，we can conclude that：The symplectic precise integration method is superior to the classic Runge-Kutta method in the numerical precision，the energy-preserving property and time-consuming；Comparing with the symplectic method，the energy can be preserved in the numerical simulation of the symplectic precise integration method.

symplectic precise integration method，Duffing equation，homogenization，symplectic

10.6052／1672-6553-2016-043

2016-07-27收到第1稿，2016-08-21收到修改稿.

*国家自然科学基金资助项目（11672233，11302169）

†通讯作者E-mail：lindu＠nwpu.edu.cn