约束多体系统动力学方程的辛算法

　　约束多体系统动力学方程的数值方法,当前是计算多体系统动力学的难点,主要是对产生的违约现象没有很好地解决,达不到数值仿真的精度要求。将通常的R2n空间的辛算法直接应用到这类动力学对应的正则形式方程,一般也要出现违约。笔者利用流形上的辛算法理论构造的一些辛算法可很好地解决违约问题。

　　1　基本理论

　　在约束多体系统动力学不考虑耗散力时,应用约束变分原理,关于完整约束g(q)=0及Hamilton量H(q,p),可得约束Hamilton正则形式方程

　　其中,q=(q1,…,qn)t∈Rn为广义坐标,p=(p1,…,pn)t∈Rn为广义动量,λ=(λ1,…,λm)t∈Rm为La-grange乘子,G(q)∈Rm×n为g(q)的Jacobian(简记为G),且秩为m,一般地假设G(q)Hpp(q,p)GT(q)可逆(它保证了系统(1)是指标为3的微分/代数方程组)。显式(隐式)一步数值积分方法,可表为,这里zi表示z(ti)的近似值。若算法Fh满足关系

　　则称Fh为辛算法。易证辛算法的合成也是辛算法。关于Runge-Kutta法,在里面可找到一类R2n空间的辛算法。现设Runge-Kutta法系数如下:

　　关于式(1)利用下列3步式辛算法^[1]

　　第1步:

　　这样构成的算法能保速度约束、位置约束不变。易证以上第1、第3步是辛算法(Euler格式的特殊情况)。所以,上面3步的合成中,若第2步也取成辛算法,则合成是辛算法,由此可构成许多约束多体系统动力学的辛算法。

　　一般的Hamilton系统在Runge-Kutta法中,记mij=biaij+bjaji-bibj, i,j=1,2,…,s,则有^[2-3]定理　若M=(mij)=0,则下面的Runge-Kutta法是辛算法

　　2　正则方程的推导^[6]

　　考查系统为带完整约束的多体系统,利用变分原理推导正则方程。

　　设系统的Lagrange函数为是广义速度,g(q)=0为系统约束,系统广义动量为则系统的Euler-Lagrange方程为

　　引入Lagrange乘子λ,又由上面两式即可得约束Hamilton系统的正则方程为

　　这样,可将以上构造的辛算法应用到以上约束正则形式的多体系统动力学方程。

　　3　算　例

　　如图1所示(N表向量误差的z-范数),有一平面三极摆,不计摆长的质量,且设摆长分别为l1=l2=l3=1 m,g=1 m/t2,m1,m2,m3分别描述每个摆的位置及集中质量,取m1=m2=m3=1 kg,初值为q1(0)=0.5,q2(0)=-,q3(0)=0,q4(0)=-2 ,q5(0)=0,q6(0)=-1-2,p1(0)=p2(0)=p3(0)=p4(0)=p5(0)=p6(0)=0,λ0=0。

　　若前面的算法的第3步中,利用二阶中点公式合成一辛算法,得仿真结果如图2、图3所示(步长0.12,时间60 s,双精度计算)。

　　若利用Baumgarte(违约校正法)法用4阶Runge-Kutta法进行仿真研究[4-6],结果则如图4、图5所示。