ⓒ 应用数学和力学编委会,ISSN 1000-0887

非线性振动分析的均向量场法

鲍四元1, 邓子辰2

(1. 苏州科技大学 土木工程学院, 江苏 苏州 215011;
2. 西北工业大学 理学院, 西安 710072)

摘要 通过构造向量形式的振动微分方程组,利用均向量场(AVF)法得到振动响应的向量差分迭代格式该离散格式能够保能量,同时具有二阶精度的特征,从而给出非线性振动问题的均向量场法介绍了均向量场法的基本步骤在建立AVF格式时,对于微分方程中若干常见的项,直接给出相应的映射项应用均向量场法研究了非线性单摆问题和Kepler(开普勒)问题,数值结果说明了该方法保能量和具有长时间求解能力的特性

均向量场法; 非线性振动; 保能量; 单摆问题; Kepler问题

引 言

动力学问题的经典算法非常丰富,如中心差分法、Runge-Kutta(龙格-库塔)法、Wilsonθ方法和Newmarkβ方法等,但这些方法本身会耗散系统的能量,其长期跟踪能力不够理想在数值计算方面,学者们正在关注和研究两大类算法:能量守恒算法和辛算法,它们都是非常优秀的算法,其中能量守恒算法能够保持系统的能量,而辛算法能够保持系统的辛结构在数值算法应尽可能多地保持原问题的本质特征指导原则下,冯康先生等首先提出辛算法的思想,他和他的研究团队取得了一系列重大成果,并且引起国内外学者的极大兴趣,产生了许多重要的后继成果[1-5]计算试验已经表明,辛算法具有优异的稳定性和长时间跟踪能力对于具体的一个守恒系统而言,用能量守恒算法积分还是用辛算法积分一般认为是比较两难的问题近年来,这一难点在一定程度上也已得到解决,由钟万勰先生及其研究团队创新性地提出了保辛-保能的算法[6]

在研究非线性动力学问题时,能保持系统能量的能量守恒算法比传统的数值算法有更大的优越性守恒系统的能量守恒格式的研究仍是现在国际前沿课题[2,7]除了经典的离散梯度法外,近来发展的平均向量场方法[8-9]就是一种构造能量守恒格式的方法(下文简称“均向量场法”)均向量场法适合于Hamilton(哈密尔顿)系统,是一种有效的离散梯度法人们对数值积分首次积分研究后发现,对于一些特殊的微分方程其数值方法的误差可以得到控制,即能够使数值方法的误差随着时间的增长而保持在很小的范围内

有一类单步法是基于B级数构造而成的,即幂级数在每一时间步长处都是向量场基本差分的和式B级数法可以保线性不变量在这一背景下,产生了保系统能量的平均向量场方法,它是一类特殊的B级数方法,可以保持所有向量场的能量

均向量场法[10-11]令人关注之处在于能够保持所有Hamilton系统的能量而自然界中的许多物理现象都可以看作Hamilton系统,其微分方程一般可表示为Hamilton微分系统国内外也有多位学者进行了相关研究,如文献[12]给出AVF方法两种局部精确和保能量的修正形式,并应用于求解Coulomb-Kepler问题;Cai等[13]给出了AVF方法在三维时域Maxwell方程中的数值分析;李昊辰等[14]研究了非线性Schrödinger(薛定谔)方程的平均向量场方法本文首先以线性微分方程为例介绍二阶精度的AVF方法,然后讨论其在两个非线性问题中的应用:单摆问题和Kepler问题针对被积函数的不同形式,给出均向量场法中具体的映射迭代格式,从而便于应用该方法

均向量场法的优势在于使用时不需要知道系统具体的能量函数,即可建立迭代格式,并能保持系统的Hamilton能量,然后通过求解代数方程得到迭代解,从而在保证精度的前提下大大方便了应用

1 均向量场法

对一般的常微分方程组

(1)

建立如下映射xnxn+1

(2)

其中h是积分步长[15]

运用内积运算和微积分中值定理等,可证明均向量场法是能够保能量的,具体参见文献[16]中的证明

式(2)是二阶精度的,也是本文所采用的格式高阶的均向量场法具有更复杂的映射形式,本文不做讨论

2 均向量场法的一般步骤

对于微分方程或微分方程组,应用均向量场法的一般步骤可分为如下3步:

1) 引入若干状态变量(往往是一些导数项),把原微分方程(组)写成向量方程组的形式;

2) 基于均向量场法写出具体的映射迭代格式;

3) 结合初值条件,迭代求解得到不同时刻的振动响应

在执行第2)步时,迭代格式右端积分项需要针对被积函数的具体形式进行化简如式(2)中,当被积函数g(x)取常见的若干函数形式时,所得积分结果分别罗列如下:

① 线性项g(u)=u,代入积分得

② 幂次项g(u)=um,积分得

③ 对正、余弦函数,积分得

④ 两个向量函数uv的乘积,积分结果为

[ξun+(1-ξ)un+1][ξvn+(1-ξ)vn+1]dξ=

在第2)步中,所得迭代格式一般是隐式的但在一些情况下可通过代数方程求解得到显式格式,需视具体问题而定

3 在线性微分系统中的应用

首先以一个线性系统为例弹簧谐振子的振动是一个线性系统,其振动方程可化为

(3)

其中x表示振动位移,y表示速度给定初始位移和速度条件x(0)=x0,y(0)=y0

由均向量场法的一般步骤,得到对应的迭代格式如下:

(4)

对式(4)进一步化简,得到一般的迭代公式:

zk+1=Azk

(5)

其中z=[x,y]′,矩阵A为Jacobi(雅可比)矩阵,形式为

其中τ为积分步长,且

对比后发现,此形式与Euler(欧拉)中点格式[17]完全相同对于Euler中点格式,文献[16-17]已经指出其在保守系统中的保辛特性及高精度由第2节中的具体映射格式知,对于线性Hamilton系统,均向量场法的迭代格式可退化为Euler-中点格式

另外,由式(4)中两式相除可得

(6)

此式恰好反映了弹簧谐振子振动过程中的能量变化情况,故该算法是保能量的

选取为计算方便,此处数值均选为量纲一图1给出该线性系统式(3)的数值模拟结果,所对比的精确解为

在数值离散方程组时时间步长的引入会导致一定的误差可选较小的时间步长τ,以减少误差实际上相对于经典算法,采用同样的时间步长时,均向量场法的误差,即使在长时间后,也是非常小的[16,18-20]

图1 一个线性微分系统的数值解与精确解的比较
Fig. 1 Comparison between the numerical solution and the exact solution for a linear differential system

图2 运动过程的相图
Fig. 2 The phase diagram of the motion

图2为τ =0.2时运动过程的相图,时间取t≤30 s由图2可知,均向量场法所得的相图始终在精确解附近的一定范围内,从而说明了均向量场法具有长时间求解的能力如果时间步长取得更小,则误差也变小

4 在单摆问题中的应用

考虑无阻尼单摆的自由振动系统[21-22],取角位移变量x,振动方程可写成

(7)

其中为单摆微幅振动固有频率的平方,l为摆长,g为重力加速度,A是初始摆角

引入新变量q=x和角速度变量方程(7)可化为

(8)

初始条件为其Hamilton函数(即系统能量)为按第2节中的方法,可写出保能量的二阶均向量场法计算公式:

(9a)

(9b)

由式(9b)/(9a),并化简得

(10)

离散时间段τ在式(10)中被消去了按单摆的力学知识,式(10)易于由动能定理得到,且其正确性与离散时间段τ的大小无关由式(10)得

(11)

式(11)中正负号根据单摆摆动方向确定,即顺时针运动取负号,而逆时针运动时为正

基于式(11)和(9b),选取若直接按固定时间步长,求解变量qp,会出现超越方程,难以求解故这里固定摆角步长dq,求出对应的不同时间间隔,从而可得动态响应(如位移、速度和加速度等)

由于单摆模型的解析解很难给出,以Runge-Kutta方法的数值解为参考比较各方法的误差图3~图6是初始摆角A为120°时的摆角-时间曲线、角速度-时间曲线,并且与ODE数值积分方法的结果(图中为实线条所示)比较,两者吻合得很好还可得出,在摆角位移随时间变化图中,接近极值点时,宜选择更小的dq,以反映真实的摆角变化

表1给出不同初始摆角时的固有频率参数T0/T的值,其中T0是线性单摆振动的周期,即由表1看出,当A=135°时,本文解(取dq=0.01 rad)与精确解的误差为0.17%;2阶摄动解的误差为6.1%;DQ法误差为0.24%而当A=150°时,本文解(取dq=0.01 rad)与精确解的误差为2.02%;2阶摄动解的误差为15.85%;DQ法误差为1.25%

图3 当dq=0.1 rad时的摆角-时间曲线 图4 当dq=0.1 rad时的角速度-时间曲线
(初摆角A=120°) (初摆角A=120°)
Fig. 3 The swing angle-time curve for Fig. 4 The angular velocity-time curve for
dq=0.1 rad (A=120°) dq=0.1 rad (A=120°)

图5 当dq=0.05 rad时的摆角-时间曲线图6 当dq=0.05 rad时的角速度-时间曲线
(初摆角A=120°) (初摆角A=120°)
Fig. 5 The swing angle-time curve for Fig. 6 The angular velocity-time curve for
dq=0.05 rad (A=120°) dq=0.05 rad (A=120°)

表1 不同初始摆角时的固有频率参数T0/T的值(dq=0.01 rad)

Table 1 Values of the natural frequency parameterT0/T with different initial swing angles(dq=0.01 rad)

5 在含耦合项的微分方程组系统中的应用

含耦合项的非线性问题有很多种,如旋转梁引起的耦合问题,又如非线性耦合振子等为简便计,这里讨论经典力学中的微分方程组模型:Kepler问题[12,19]其对应的微分方程组为

(12a)

(12b)

其中q1,q2的物理意义是质点位置的横、纵坐标式(12)可写成如下方程组:

(13)

其Hamilton函数是守恒量,而且,角动量Θ=q1p2-q2p1也是守恒量

基于二阶均向量场法的步骤,得Kepler问题中式(13)的等效数值格式如下:

(14a)

(14b)

(14c)

(14d)

分别运算式(14a)/(14b)和式(14c)/(14d),可消去τ,得

(15a)

(15b)

由式(15a)+(15b),化简得

(16)

式(16)即对应于原系统的Hamilton量(能量)守恒这也间接地验证了均向量场法格式是保能量的

图7给出时间步长取不同值(0.03 s和0.07 s)时的轨道图形 当dt取0.07 s时,所得结果在图7(a)上显示有一定的进动,这对于数值积分而言误差是难免的;而小步长0.03 s时,在长时间内几乎没有出现进动图8是积分时间为18 s时,本文解与精确解的比较由图8可知,取大步长0.07 s时,轨道图的初始误差较小,然后逐步增大但步长取0.02 s时,本文解与真实椭圆轨道解基本重合,长时间内的误差也很小

图9给出dt=0.07 s时的角动量误差随时间的变化, 其中初始角动量为Θ0=0.8, 此处的误差定义为(Θ-Θ0)/Θ0,不同时刻的误差均在10-2以内图10给出了系统能量随时间的变化, 其中初始Hamilton函数H0=0.5,能量误差定义为(H-H0)/H0,不同时刻的误差均在10-5以内

(a) dt=0.07 s(b) dt=0.03 s
图7 不同dt时的轨道图
Fig. 7 The orbit diagrams with different dt values

(a) dt=0.07 s(b) dt=0.02 s
图8 本文解与轨道真实解[24]间的比较
Fig. 8 Comparison of orbit diagrams between the present solution and the exact solution[24]

图9 当dt=0.07 s时的角动量误差图10 当dt=0.07 s时Hamilton守恒量误差
随时间的变化 随时间的变化
Fig. 9 The angular momentum error varying Fig. 10 The Hamiltonian conserved quantity error
with time, dt=0.07 s varying with time, dt=0.07 s

6 结 论

1) 采用保能量的均向量场法研究非线性振动方程,需先转为状态空间的一阶微分方程组,利用映射格式建立迭代方程组,该方法能够保持Hamilton系统的能量

2) 本文给出了常见函数在均向量场格式中积分后的转化形式,便于直接应用该法

3) 对于非线性程度较高的单摆振动,甚至含耦合项的微分振动系统,算例结果说明本文方法在积分步长相对较大时精度较二阶摄动法高,与微分求积法精度接近如果减少积分步长,精度会更好

4) 在解决线性振动系统时,均向量场法可退化为经典的Euler-中点法

参考文献

[1] 冯康, 秦孟兆. 哈密尔顿系统的辛几何算法[M]. 杭州: 浙江科学技术出版社, 2003.(FENG Kang, QIN Mengzhao.Symplectic Geometric Algorithm for Hamilton System[M]. Hangzhou: Zhejiang Science and Technolog Press. 2003.(in Chinese))

[2] 秦孟兆, 王雨顺. 偏微分方程中的保结构算法[M]. 杭州: 浙江科技出版社, 2010.(QIN Mengzhao, WANG Yushun.Structure-Preserving Algorithms for Partial Differential Equation[M]. Hangzhou: Zhejiang Science and Technolog Press, 2010.(in Chinese))

[3] FENG K. Difference schemes for Hamiltonian formalism and symplectic geometry[J].Journal of Computational Mathematics, 1986,4(3): 279-289.

[4] 胡伟鹏, 邓子辰. 无限维动力学系统的保结构分析方法[M]. 西安: 西北工业大学出版社, 2015.(HU Weipeng, DENG Zichen.The Infinite Dimensional Dynamical System Structure Analysis Method[M]. Xi’an: Northwestern Polytechnical University Press, 2015.(in Chinese))

[5] HU W P, DENG Z C, ZHANG Y. Multi-symplectic method for peakon-antipeakon collision of quasi-Degasperis-Procesi equation[J].Computer Physics Communications, 2014,185(7): 2020-2028.

[6] 高强, 钟万勰. Hamilton系统的保辛-守恒积分算法[J]. 动力学与控制学报, 2009,7(3): 193-199.(GAO Qiang, ZHONG Wanxie. The symplectic and energy preserving method for the integration of Hamilton system[J].Journal of Dynamics and Control, 2009,7(3): 193-199.(in Chinese))

[7] BRUGNANO L, IAVERNARO F, TRGIANTE D. A two-step, fourth-order method with energy preserving properties[J].Computer Physics Communications, 2012,183(9): 1860-1868.

[8] 陈璐, 王雨顺. 保结构算法的相位误差分析及其修正[J]. 计算数学, 2014,36(3): 271-290.(CHEN Lu, WANG Yushun. Phase error analysis and correction of structure preserving algorithms[J].Mathematica Numerica Sinica, 2014,36(3): 271-290.(in Chinese))

[9] 叶霄霄. 基于平均向量场方法的暂态稳定计算[D]. 硕士学位论文. 宜昌: 三峡大学, 2015.(YE Xiaoxiao. Transient stability calculation based on the average vector field method[D]. Master Thesis. Yichang: China Three Gorges University, 2015.(in Chinese))

[10] QUISPEL G R W, MCLACHLAN D I. A new class of energy-preserving numerical integration methods[J].Journal of Physics A:Mathematical and Theoretical, 2008,41(4): 045206. DOI: 10.1088/1751-8113/41/4/045206.

[11] CELLEDONI E, MCLACHLAND I, OWREN B, et al. Energy-preserving integrators and the structure of B-series[J].Foundations of Computational Mathematics, 2010,10(6): 673-693.

[12] J L. Improving the accuracy of the AVF method[J].Journal of Computational and Applied Mathematics, 2014,259: 233-243.

[13] CAI J X, WANG Y S, GONG Y Z. Numerical analysis of AVF methods for three-dimensional time-domain Maxwell’s equations[J].Journal of Scientific Computing, 2016,66(1): 141-176.

[14] 李昊辰, 孙建强, 骆思宇. 非线性薛定谔方程的平均向量场方法[J]. 计算数学, 2013,35(1): 60-66.(LI Haochen, SUN Jianqiang, LUO Siyu. An averaged vector field method for the nonlinear Schrödinger equation[J].Mathematica Numerica Sinica, 2013,35(1): 60-66.(in Chinese))

[15] HAIRER E, LUBICH C, WANNER G.Geometric Numerical Integration:Structure-Preserving Algorithms for Ordinary Differential Equations[M]. Berlin: Springer, 2006.

[16] 陈璐. 保结构算法的相位误差分析及其修正[D]. 硕士学位论文. 南京: 南京师范大学, 2014.(CHEN Lu. Phase error analysis and correction of structure preserving algorithms[D]. Master Thesis. Nanjing: Nanjing Normal University, 2014.(in Chinese))

[17] 邢誉峰, 杨蓉. 动力学平衡方程的中点辛差分求解格式[J]. 力学学报, 2007,39(1): 100-105.(XING Yufeng, YANG Rong. Application of Euler midpoint symplectic integration method for the solution of dynamic equilibrium equations[J].Chinese Journal of Theoretical and Applied Mechanics, 2007,39(1): 100-105.(in Chinese))

[18] 刘晓梅, 周钢, 王永泓, 等. 辛算法的纠飘研究[J]. 北京航空航天大学学报, 2013,39(1): 22-26.(LIU Xiaomei, ZHOU Gang, WANG Yonghong, et al. Rectifying drifts of symplectic algorithm[J].Journal of Beijing University of Aeronautics and Astronautics, 2013,39(1): 22-26.(in Chinese))

[19] 邢誉峰, 杨蓉. 单步辛算法的相位误差分析及修正[J]. 力学学报, 2007,39(5): 668-671.(XING Yufeng, YANG Rong. Phase errors and their correction in symplectic implicit single-step algorithm[J].Chinese Journal of Theoretical and Applied Mechanics, 2007,39(5): 668-671.(in Chinese))

[20] 秦于越, 邓子辰, 胡伟鹏. 谐振子的辛欧拉分析方法[J]. 动力学与控制学报, 2014,12(1): 9-12.(QIN Yuyue, DENG Zichen, HU Weipeng. Symplectic Euler method for harmonic oscillator[J].Journal of Dynamics and Control, 2014,12(1): 9-12.(in Chinese))

[21] 李鹏松, 孙维鹏, 吴柏生. 单摆大振幅振动的解析逼近解[J]. 振动与冲击, 2008,27(2): 72-74.(LI Pengsong, SUN Weipeng, WU Baisheng. Analytical approximate solutions to large amplitude oscillation of a simple pendulum[J].Journal of Vibration and Shock, 2008,27(2): 72-74.(in Chinese))

[22] 吕中荣, 刘济科. 摆的振动分析[J]. 暨南大学学报(自然科学版), 1999,20(1): 42-45.(LÜ Zhongrong, LIU Jike. Vibration analysis of a pendulum[J].Journal of Jinan University(Natural Science), 1999,20(1): 42-45.

[23] 周凯红, 王元勋, 李春植. 微分求积法在单摆非线性振动分析中的应用[J]. 力学与实践, 2003,25(3): 50-52.(ZHOU Kaihong, WANG Yuanxun, LI Chunzhi. The application of differential quadrature method in nonlinear vibration analysis of simple pendulum[J].Mechanics in Engineering, 2003,25(3): 50-52.(in Chinese))

[24] 李文博, 赵定柏. 开普勒问题的一种简单处理[J]. 大学物理, 2000,19(1): 45-47.(LI Wenbo, ZHAO Dingbai. A simple treatment of the Kepler problem[J].College Physics, 2000,19(1): 45-47.(in Chinese))

An Average Vector Field Method for Nonlinear Vibration Analysis

BAO Siyuan1, DENG Zichen2

(1.School of Civil Engineering,Suzhou University of Science and Technology,
Suzhou,Jiangsu 215011,P.R.China
2.School of Natural and Applied Sciences,
Northwestern Polytechnical University,Xian 710072,P.R.China)

(Contributed by DENG Zichen, M. AMM Editorial Board)

Abstract: Through construction of differential equations in the vector form, the differential iteration form of the vibration response was obtained according to the average vector field (AVF) method. This discrete form is energy-preserving for the Hamiltonian system, and has the characteristics of 2nd-order accuracy. The detailed steps of the AVF method were given. To establish the AVF scheme, the mapping forms were deduced directly for several common items in the differential equations. The pendulum problem and the Kepler problem were studied with the AVF method. The numerical results demonstrate the advantages of the AVF method in solving nonlinear vibration problems, i.e. the conservation of energy and the long-term solution stability.

Key words: average vector field method; nonlinear vibration; energy-preserving; pendulum problem; Kepler problem

Foundation item: The National Natural Science Foundation of China(11202146)

(我刊编委邓子辰来稿)

中图分类号 O322; O326

文献标志码:A

DOI: 10.21656/1000-0887.390178

文章编号1000-0887(2019)01-0047-11

收稿日期 2018-06-26;

修订日期2018-11-07

基金项目 国家自然科学基金(11202146)

作者简介 鲍四元(1980—),男,副教授(E-mail: bsiyuan@126.com);邓子辰(1964—),男,教授,博士生导师(通讯作者. E-mail: dweifen@nwpu.edu.cn).

引用本文/Cite thispaper:鲍四元, 邓子辰. 非线性振动分析的均向量场法[J]. 应用数学和力学, 2019,40(1): 47-57.BAO Siyuan, DENG Zichen. An average vector field method for nonlinear vibration analysis[J].Applied Mathematics and Mechanics, 2019,40(1): 47-57.