基于球摆模型的离散变分积分子算法研究*

白龙 戈新生

(北京信息科技大学机电工程学院，北京 1 00192)

基于球摆模型的离散变分积分子算法研究*

白龙†戈新生

(北京信息科技大学机电工程学院，北京 1 00192)

在动力学系统长时间的仿真计算中，力学系统固有的结构将影响到计算精度及稳定性．离散变分积分子能够保持力学系统的能量，动量及辛结构的守恒．结合离散变分原理，通过对系统的拉格朗日函数进行离散化以及求变分和积分的过程，可以得到力学系统的离散变分积分子算法．该算法是一种递归算法，给定初始条件便可得到系统的动力学参数的时间历程．使用该原理可以构造具有完整约束的拉格朗日系统的辛－动量积分子方法．与连续算法相比，离散变分积分子算法能够直接在离散拉格朗日函数的基础上得到姿态与角速度的递推公式，而不需要复杂的迭代计算．本文研究是基于第一类拉格朗日函数的离散变分积分子算法．球摆模型是一个具有完整约束的拉格朗日系统．仿真结果表明，系统的能量值在长时间的仿真中得到保持，且计算的精度与步长的数量级呈现二次方的关系，系统角速度和姿态的仿真结果都符合球摆的运动规律．

离散变分原理， 力学积分子， 拉格朗日函数， 能量守恒

引言

在动力学系统的实际计算中，能量守恒与动量守恒与算法能否保持系统的结构具有很大关系．通过拉格朗日函数得到系统的运动微分方程的解法在长时间的仿真中不具有保持系统结构的性质，而力学积分子能够保持力学系统结构．因此从理论上来讲，它能够保持系统的能量，动量及辛结构的守恒．通过离散变分原理，得到离散拉格朗日方程，进而得到力学系统的离散变分积分子．与传统算法比较，该方法既能保证计算的效率，又具有较好的长时间仿真的性质．

离散变分积分算法经过近半个世纪的发展，已经成为一个较完整的算法体系．其中Verlet算法［1］可以用来处理无约束的力学系统，Shake算法［2］是Verlet算法的延伸，用来处理带有完整约束的力学系统．为处理完整约束系统因离散而带来的速度约束，Rattle算法［3］作为Shake算法的延伸较好的处理了该问题．

近年来，基于李群的离散变分积分子算法得到广泛研究［4］，并且在3D刚体摆的动力学计算中取得了良好的效果．该理论摆脱了传统离散变分积分子算法的约束理论，利用李群理论保证了系统的结构，使系统的能量和动量在长时间的仿真中都保持守恒．

本文研究的是带有约束的拉格朗日系统的保能量离散变分积分子方法．该方法通过离散变分积分子方法，得到一个离散的拉格朗日函数，使得约束在每个时间步长都得到满足．在长时间的仿真中，离散拉格朗日方程保证了辛结构和动量的守恒，进而保证了能量的守恒．球摆是3D摆的一种简化模型，在离散变分积分子算法下，其角速度、姿态是否符合其运动规律将为离散变分积分子算法应用到3D摆的研究中提供理论保证．仿真结果表明，在离散变分积分子算法下，系统的能量维持恒定，并随着时间步长的减小，计算的精度也逐步提高．

1 离散变分原理

给定一个结构空间Q，离散拉格朗日函数是一个映射L:Q×→．

对于固定的正整数N，对L(qk+1，qk)求和运算是映射S:QN+1→ ，定义为

其中，qk∈Q，k∈ 是离散时间．对 L(qk+1，qk)求和是对L(q，q)积分的一个离散近似．通过对S在q1，…，qN－1求极限，得到离散形式的欧拉—拉格朗日方程

其中k∈{1，…，N－1}，Φ:Q×Q→Q×Q定义为 Φ(qk，qk－1)=(qk+1，qk)．如果D2L 是可逆的，那么方程(3)定义了离散映射Φ，使得系统在离散时间下向前演化．定义映射

通过使用T*Q上的规范二阶形式，定义Q×Q上的二阶量ω

选择Q上的坐标qi和T*Q上的正则坐标(qi，pi)，该坐标满足 ΩCAN=dqi∧dpi以及 ΘCAN=pidqi，那么离散欧拉－拉格朗日方程是:

根据映射 Φ(qk+1，qk)=(qk+1，qk)得到

对方程(5)计算得到:

由于∧=0，则

由式(7)，对ω求Φ运算

由式(9)可知:Φ－1ω=ω，即 Φ 能够保持系统的辛流形．

令离散拉格朗日函数在Q上的李群G的对角线作用下恒定，令ξ∈g，其中g是G的李代数．则

对式(11)求微分，令s=0得到:

其中，ξQ是无穷小生成元．考虑对L(qk+1，qk)的求和运算，(0 ＜k＜N)，通过qk+1(s)=exp(sξ)qk+1将qk+1变s∈ 的形式．根据qk+1(0)对S求极限得到:

由方程(7)，(12)还可以写为以下形式:

从方程(14)减去方程(12)得到:

如果动量映射J:Q×Q→g*定义为:

则方程(16)说明动量映射由Φ:Q×Q→Q×Q保持．

2 离散拉格朗日方程

球摆模型如图1所示，质心在惯性坐标系中的投影可以表示为

图1 球摆模型Fig．1 The model of spherical pendulum

球摆的动能为

势能可以写为

由拉格朗日函数L=T－V可得

约束条件为

定义表达式:

将q和q的离散表达式分别代入连续拉格朗日函数中，得到离散拉格朗日函数:

则连续拉格朗日函数(21)的离散形式为:

由离散勒让德变换可知:

分别在qk，qk+1对约束条件g(q)=0进行离散化可以得到

使用拉格朗日乘子，将约束条件添加到离散拉格朗日函数后，求变分可以得到

在考虑因离散化而对拉格朗日系统产生不显含速度约束的情况时，离散系统在不显含速度的作用下，会使系统产生耗散现象，导致仿真结果不具有一般性，因此，要解决此问题，就要在离散动量递推公式后再添加一个约束项，对系统的速度进行约束．

由 Rattle算法［3］的定义

将式(31)代入式(30)得到qk+1的递推公式．再将得到的递推公式代入约束条件g(qk+1=0)可求出式(31)中常量γk．

由式(26)可以得到系统的动量的递推公式表示为

对约束条件g(qk+1)=0求导得到包含速度的约束条件

将式(30)、(32)代入式(33)便可求出未知量λk．

3 仿真实验

球摆模型是3D摆的一种简化模型，设摆的重量为 1 kg，g=9．8N/kg，摆长 l =1m，选取初始位移坐标q=()以及初始速度q=(0，0，0)，时间步长为 0．001s．仿真结果如图所示．

图2和图3表明，球摆系统随时间呈现周期性运动，由于q和q表征的是系统在惯性坐标系下的运动，根据式(17)和(18)得出θ和φ以及θ·和 φ随时间的变化图，如图4、图5所示．

图2 球摆坐标的变化规律Fig．2 The coordinate variation of the spherical pendulum

图3 球摆速度的变化规律Fig．3 The velocity variation of the spherical pendulum

图4 球摆角度的变化规律Fig．4 The angle variation of the spherical pendulum

图5 球摆角速度的变化规律Fig．5 The angle velocity variation of the spherical pendulum

图4表明，球摆的θ角度呈现周期性变化;φ角保持恒定，符合球摆的运动规律．

图5表明，球摆的角速度也呈现周期性变化，且 φ在较高的精度下保持恒定．

图6 球摆动能势能及总能量图Fig．6 The kinetic potential and total energy of spherical pendulum

由图6可知，球摆系统的能量值呈现周期性变化，符合球摆运动时能量的变化规律，而且系统的总能量在宏观上维持恒定，未出现随时间增加而逐渐偏离理论值的现象．

图7 球摆的总能量比较Fig．7 the comparison of the total energy of the spherical pendulum

由图7可知，与四阶龙格库塔算法相比，离散变分积分子算法能够保持系统的总能量维持在一恒定的范围波动，不会出现随仿真时间的增加而逐渐发散的现象．

表1为在不同数量级的步长下，系统的能量的波动范围及误差随时间步长的变化．通过对比时间步长与误差可以得出，误差的数量级与时间步长的数量级呈现二次方的关系．即误差随步长的减小而逐渐减小．

表1 系统总能量变化Table 1 The fluctuate of the energy

4 总结与展望

离散变分积分子算法作为一种保结构的算法，在力学系统长时间的仿真中有着明显的优势．本文表明，在简单非线性系统的计算中，传统的离散变分积分子算法取得了长时间保能量的良好的效果，如何将其应用到更加复杂的诸如3D刚体摆等非线性系统的计算中将是下一步研究的重点．

传统离散变分积分子算法是基于类似于带有约束条件的第一类拉格朗日方程的形式来达到保能量的要求，3D刚体摆是具有三个转动自由度的完整系统，欧拉四元数的建模方法为将离散变分积分子算法应用到3D刚体摆中提供了可能，该方法通过四个变量实现对3D刚体摆的描述，而四个变量满足的约束条件可以充当3D刚体摆在应用离散变分积分子算法时所需的约束条件．

除传统离散变分积分子算法以外，基于李群的离散变分积分子算法通过系统运动学方程中的反对称结构来保持系统的辛流形，从而达到保能量的要求．该算法因其不需要单独寻找约束条件而更加简洁有效．国内外相关研究也证明，基于李群的离散变分积分子算法也可应用到3D刚体摆的动力学计算中．

1 Anderson H C．Rattle:A velocity version of the shake algorithm for molecular dynamics calculations．Journal of Computational Physics，1983，52:24 ～34

2 Wendlandt J M，Marsden J E．Mechanical integrators derived from a discrete variational principle．Physica D:Nonlinear Phenomena，1997:223～246

3 Bacconi，B．Spacecraft attitude dynamics and control［PhD Thesis］．Michigan:Department of Aerospace Engineering University of Michigan．2006:25 ～33

4 Wendlandt J M．Control and simulation of multibody systems［PhD Thesis］．Berkeley:University of California．1997:18～40

5 Shen J L．Dynamics and control of a 3 D pendulum．Proceedings of the IEEE Conference on Decision and Control，2004 43rd IEEE Conference on Decision and Control(CDC)，2004:323 ～328

6 Chaturvedi N A．Nonlinear dynamics of the 3D pendulum．Journal of Nonlinear Science，2011:3 ～32

7 丁洁玉，潘振宽．非完整约束多体系统时间离散变分积分法．动力学与控制学报，2011，9(4):289～292(Ding J Y，Pan Z K．Time discrete variational integrator for multibody dynamic systems with nonholonomic constraints．Journal of Dynamics and Control，2011，9(4):289 ～ 2 92)(in Chinese)

8 Marsden J E，West M．Discrete mechanics and variational integrators．Acta Numerica，2001:357 ～514

9 Lee，Taeyoung．A lie group variational integrator for the attitude dynamics of a rigid body with applications to the 3D pendulum．Proceeding of the IEEE International Conference on Control Applications，2005:962 ～ 9 67

10 Lee T Y．Lie group variational integrators for the full body problems in orbital mechanics．Celestial Mech Dyn Astr，2007，98:121 ～144

11 丁希仑，刘颖．用李群李代数分析具有空间柔性变形杆件的机器人动力学．机械工程学报，2007，43(12):184～189(Ding X L，Liu Y．Dynamic analysis of robot with spatial compliant links using lie group and lie algebra．Chinese Journal of Mechanical Engineering，2007，43(12):184～189(in Chinese))

12 高强，钟万勰．非完整约束动力系统的离散积分方法．动力学与控制学报，2012，10(3):193～198(Gao Q，Zhong W X．Numerical algorithms for dynamic system with non －holonomic constraints．Journal of Dynamics and Control，2012，10(3):193 ～198(in Chinese))

13 Mclachlan R，Perlmutter M．Integrators for nonholonomic mechanical systems．Journal of Nonlinear Science，2006，16:283～332

14 钟万勰，高强．约束动力系统的分析结构力学积分．动力学与控制学报，2006，4(3):193～200(Zhong W X，Gao Q．integration of constrained dynamical system via analytical structual mechanicals．Journal of Dynamics and Control，2006，4(3):193 ～200(in Chinese))

*The project supported by the National Natural Science Foundation of China(11072038)the B category of State key development program of the Natural Science Foundation of Beijing(KZ201110772039)

† Corresponding author E－mail:bailong0316jn@126．com

THE DISCRETE VARIATIONAL INTEGRATORS METHOD OF THE SPHERICAL PENDULUM*

Bai Long†Ge Xinsheng
Beijing Information Science＆Technology University，Beijing100192，China

The mechanical system＇s intrinsic structure may influence the long time computation＇s accurate and stability．The discrete variation integrators can conserve the energy momentum and the symplectic structure of the system．Combined with the discrete variation principle，the discrete variation integrator method can be obtained through the process of discretization varaiton and integration．This is a recursive algorithm that the time history of the parameters only need the initial condition．According this theory，a sympletic－momentum integrator can be formulated for the holonomic constraint Lagrange system．This method can get the recursive formula of the attitude and the angle velocity direct form the discrete Lagrange function and don＇t need complicated iterative computation．The discrete variation integrator method explored in this paper is based on the first Lagrange function．The spherical pendulum is a Lagrange system with holonomic constraints．The simulation result states that the energy be conserved in a long time simulation，and the accuracy of the computation presents a quadratic relation with the time step．The angle velocity and the attitude also present different character under two different algorithm．

discrete variation principle， mechanical integrator， Lagrange function， energy conservation

21 August 2012，

29 May 2013．

10．6052/1672－6553－2013－077

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

*国家自然科学基金资助项目(11072038)、北京市自然科学基金重点项目B类(KZ201110772039)

E－mail:bailong0316jn@126．com