谐振子的辛欧拉分析方法*

秦于越 邓子辰 胡伟鹏

(1．西北工业大学力学与土木建筑学院，西安 710072)(2．大连理工大学工业装备结构分析国家重点实验室，大连 116023)

谐振子的辛欧拉分析方法*

秦于越1†邓子辰1，2胡伟鹏1

(1．西北工业大学力学与土木建筑学院，西安 710072)(2．大连理工大学工业装备结构分析国家重点实验室，大连 116023)

针对理想简谐振子力学模型，研究了其守恒律，并利用辛欧拉格式分析简谐振子振动过程．首先给出了谐振子系统的平方守恒律、周期守恒律和相差守恒律．构造了谐振子的普通欧拉格式和辛欧拉格式，研究了两种格式下三种守恒律各自的保持情况．模拟结果显示:辛欧拉格式能够精确保持时域守恒律(平方守恒律)，但无法保持频域守恒律(周期守恒律和相差守恒律)．如要克服辛欧拉格式的不足，需按邢誉峰教授提出的方法进行校正．

哈密尔顿， 保结构算法， 辛欧拉， 简谐振子， 守恒律

引言

简谐振子是固体力学领域最简单的理想力学模型之一．由于忽略系统的阻尼等耗散因素，简谐振子的振动过程是一个典型的保守系统，振动能量应该是守恒的，这一点在理论上早已得到了学术界的公认，但是从数值分析中研究这类问题具有更广阔的工程应用背景．

1984年，保结构几何积分算法创始人在双微国际会议上作了题为“On difference schemes and symplectic geometry”的大会报告［1］，首次为物理力学领域的数值算法提出了在数值求解过程中应该充分关注物理力学系统的固有几何性质的保持，得到了国际同行的广泛认同，冯康先生也因此获得了1997年自然科学一等奖．此后，辛算法在诸多领域得到了应用:钟万勰院士将辛算法引入应用力学领域的研究，并创造性地提出了精细积分算法用于时间历程的积分运算［2］;秦孟兆等人继承了冯康先生辛算法思想，提出了多步辛算法的构造方法，并将其应用于波动问题的数值分析过程，得到了一些有意义的结果［3，4］;Hairer教授等人全面总结辛算法研究成果，系统提出了针对常微分方程的保结构几何积分算法理论体系［5］;Budd 和 Piggott在 Hamilton框架下，系统阐述了保结构算法的基本思想，为后续保结构算法的发展指明了方向［6］;正是受保结构算法思想的启发，Bridges教授等人将针对有限维Hamilton动力学系统的辛算法推广至无限维，创立了多辛算法［7］;邓子辰教授等人针对非保守系统给出了广义多辛算法［8－11］．

本文针对简谐振子模型，研究了简谐振子的固有几何性质及其对应的守恒律，构造辛欧拉格式模拟简谐振子的振动过程，通过与常规欧拉格式的模拟结果进行对比，充分展现辛欧拉格式的保结构特点和长时间数值稳定性特点．

1 简谐振子的守恒律

考虑以下简谐振子谐振问题

这是一个极其简单的数学物理问题，其精确解可以很容易地得到，其中一个最为简单的精确解为:

这个精确解很显然是周期有界的，结合代数理论，得到其满足的守恒律如下:

①平方守恒律:

这一守恒律表明谐振子的振动应该始终都在一个闭合的圆上．

②周期守恒律:

其中T(u)和T(v)分别表示两个广义位移分量的振动周期．这一守恒律表明谐振子的两个广义位移分量的振动周期均是2π;

③相差守恒律:

需要说明的是以上守恒律是不依赖时间的，也就是说，按照保结构思想，经过无限长的时间后，以上三个守恒律仍然是应该得到精确满足的．

2 谐振问题的辛欧拉离散格式

令时间步长为h，依照差分理论，谐振问题(1)的普通欧拉差分离散格式为:

与严格保能量的欧拉中点辛差分［12］格式相比，格式(6)是一个显式格式，计算速度会更快．从这一差分格式中很容易就可以得到第n+1步的结果与第n步的结果之间的代数关系为:

这一代数关系不仅不依赖于时间，而且不依赖于初始条件和边界条件．代数关系(7)从根本上与平方守恒律(3)相背．同时，无论时间步长多么小，由普通欧拉格式(6)得到的谐振子的振动将是无界的，即，随着时间的演化，振动幅值趋近于无穷大，这从根本上破坏了原系统的周期性，更谈不上数值解能满足周期守恒律(4)和相差守恒律(5)了．

注意到系统(1)是一个典型的分离形式的Hamilton系统，其Hamilton数为:

采用辛欧拉离散方法［5］，离散系统(1)，得到离散辛格式:

其迭代雅克比矩阵为:

容易验证该雅克比矩阵是辛的．

3 数值试验

假定谐振子的初始条件为:

图1 谐振子两个广义位移分量变化情况Fig．1 The evolutions of functions and

取时间步长h=0．05s，分别采用普通欧拉格式(6)和辛欧拉格式(9)模拟谐振子的振动过程，得到两个广义位移分量结果如图1，同时得到三个守恒律的保持情况如图2和图3(后两种守恒律的保持情况只能给出辛欧拉格式的结果，因为普通欧拉格式得到的结果是非周期的)．

图2给出了两种格式得到的数值结果在平方守恒律方面的结果对比，从对比结果不难发现:普通欧拉格式得到的结果满足方程(7)，即u2+v2的数值线性增长;而辛欧拉格式的u2+v2数值始终保持为1．

图3 周期守恒律和相差守恒律结果Fig．3 The conservation law of period and phase difference

从周期守恒律和相位差守恒律结果可以看出:随着模拟的进行，谐振子两个广义位移分量的周期缓慢线性增长，同时两个广义位移之间的相位差也缓慢线性增长，这说明辛算法不能精确保持系统的频域特性．这一点与邢誉峰教授关于辛算法的相位误差的研究结果［13］吻合，如果要精确保持系统的频域特性，需要按照邢誉峰教授建议的方法进行校正．

4 结论

本文针对谐振子理想模型，分析其多种守恒律，并从理论和数值模拟结果两个方面分别比较了谐振子模型普通欧拉格式和辛欧拉格式在保持守恒律方面的差异．对比结果表明:辛欧拉格式能够较好的长时间模拟谐振子的振动过程，同时长时间精确保持谐振子时域内的平方守恒律，但是不能精确保持频域内的周期守恒律和相位差守恒律，若要求其能够精确保持频域内的相关守恒律，还需要对辛欧拉格式进行校正．

1 Feng K．On difference schemes and symplectic geometry，Proceeding of the 1984 Beijing symposium on D．D．，Beijing:Science Press，1984，42～58

2 Zhong W X．Some developments of computational solid Mechanics in China．Computers＆Structures，1988，30(4):783～788

3 Qin M Z，Zhang M Q．Multi－stage symplectic schemes of two kinds of Hamiltonian systems for wave equations．Computers＆Mathematics with Applications，1990，19(10):51～62

4 Zhao P F，Qin M Z．Multisymplectic geometry and multisymplectic Preissmann scheme for the KdV equation．Journal of Physics a－Mathematical and General，2000，33(18):3613～3626

5 Hairer E，Lubich C，Wanner G．Geometric numerical integration:structure preserving algorithms for Ordinary Differential Equations．Berlin:Springer－Verlag，2002

6 Budd C J，Piggott M D，Geometric integration and its applications．Amsterdam:Handbook of Numerical Analysis，2003:35～139

7 Bridges T J．Multi－symplectic structures and wave propagation．Mathematical proceedings of the Cambridge Philosophical Society，1997，121(1):147～190

8 Hu W P，Deng Z C，Han S M，Zhang W R．Generalized multi－symplectic integrators for a class of Hamiltonian nonlinear wave PDEs．Journal of Computational Physics，2013，235:394～406

9 胡伟鹏，邓子辰．大阻尼杆振动的广义多辛算法．动力学与控制学报，2013，11(1):1～4(Hu W P，Deng Z C．Generalized multi－symplectic method for vibration of big damping bar．Journal of Dynamics and Control，2013，11(1):1～4(in Chinese))

10 Hu W P，Deng Z C，Wang B，Ouyang H J．Chaos in an embedded single－walled carbon nanotube．Nonlinear Dynamics，2013，72(1 －2):389 ～398

11 Hu W P，Deng Z C，Ouyang H J．Generalized Multisymplectic Method for Dynamic Responses of Continuous Beam under Moving Load．International Journal of Applied Mechanics，2013，5(3):1350033

12 邢誉峰，杨蓉．动力学平衡方程的Euler中点辛差分求解格式．力学学报，2007，39(1):100～105(Xing Y F，Yang R．Application of Euler midpoint symplectic integration method for the solution of dynamic equilibrium equations．Acta Mechanica Sinica，2007，39(5):100～105(in Chinese))

13 邢誉峰，杨蓉．单步辛算法的相位误差分析及修正．力学学报，2007，39(5):668～671(Xing Y F，Yang R．Phase errors and their correction in symplectic implicit sin－gle－step algorithm．Acta Mechanica Sinica，2007，39(5): 668 ～671(in Chinese))

*The project supported by the National Natural Science Foundation of China(11172239，11002115，11372253)，the Doctoral Program Foundation of Education Ministry of China(20126102110023)and the Open Foundation of State Key Laboratory of Structural Analysis of Industrial Equipment(GZ0802)

† Corresponding author E－mail:769482448@qq．com

SYMPLECTIC EULER METHOD FOR HARMONIC OSCILLATOR*

Qin Yuyue1†Deng Zichen1，2Hu Weipeng1
(1Department of Engineering Mechanics，Northwestern Polytechnical University，Xi’an710072，China)(2State Key Laboratory of Structural Analysis of Industrial Equipment，Dalian University of Technology，Dalian116023，China)

Focusing on the conservation properties，the symplectic Euler scheme of the harmonic oscillator was constructed to analyze its vibration properties．Firstly，three conservation laws，including the square conservation law，the period conservation law and the phase difference conservation law，were presented for the harmonic oscillator．And then，the common Euler scheme and the symplectic Euler scheme were constructed to study the above three conservation laws．The numerical results imply that the symplectic Euler scheme can preserve the conservation law in time domain(the square conservation law)exactly，but can’t preserve the conservation laws in phase domain(the period conservation law and the phase difference conservation law)，which is the shortcoming of the symplectic method but can be overcome by the modification method presented by Prof．Xing．

Hamiltonians， structure － preserving method， symplectic Euler， harmonic oscillator， conservation law

21 January 2013，

29 August 2013．

10．6052/1672－6553－2013－094

2013－01－21 收到第 1 稿，2013－08－29 收到修改稿．

*国家自然科学基金(11172239，11002115，11372253)，博士点基金(20126102110023)和大连理工大学工业装备结构分析国家重点实验室开放基金(GZ0802)

E－mail:769482448@qq．com