（天体物理专业论文）哈密顿摄动分解体系的四阶力梯度辛算法.pdf

摘要 摘要 描述哈密顿系统的混沌运动需要依赖于可靠的数值方法和混沌识别方法。 本文主要工作在于数值方法的开发与利用。 我们在c h i n 等提出的四阶力梯度辛算法基础上构造了两类显式四阶力梯度 辛算法；进一步将这些力梯度辛算法推广应用到由w i s d o m 和h o l m a n 对太阳系 n 体问题所给出的哈密顿函数摄动分解情形；最后利用推广形式力梯度辛算法研 究行星磁气圈内的带电粒子的混沌运动。分别简述如下。 首先，对于分解为二次动能r 和势能y 两部分的哈密顿系统，在势能y 对应 的l i e 算子中加入力梯度算子通过对称组合我们构造了一组显式四阶力梯度辛 算法，其中包括c h i n 等提出的力梯度辛算法。这些算法能够推广用来求解j a c o b i 坐标下所分解成的k e p l e r 部分日。和摄动部分h ，的n 体引力哈密顿系统。数值 结果表明，在r + v 型哈密顿分解下，每个梯度算法的精度大大优越于f o r e s t r u t h 的非力梯度四阶辛方法；但是对于h 。+ 日。型分解情形，就平经度和相对位置精 度而言每个梯度算法与f o r e s t r u t h 的四阶辛方法几乎等效。同时，无论是在 丁+ v 型分解还是在风+ h ，型分解中，这些梯度算法在数值性能上没有明显差 异。应当着重指出的是后种分解与前种分解相比每个梯度算法能大大提高数值 精度。由于这种推广具有快速和高精度的优点，所以值得推荐用来模拟n 体问 题的各种轨道运动。 其次，利用新提出的四阶力梯度辛算法并结合日。+ 日型摄动分解情形我们 数值研究了行星磁气圈内的带电粒子的混沌运动。该物理模型可以简化为一个 摄动二体哈密顿问题。带电粒子在赤道平面上的几种相图形状取决于电荷质量 比和z 方向角动量这两个动力学参数，当然离开赤道平面的运动除这两个动力学 参数外还受能量的影响。发现当增大能量或电荷质量比的绝对值时混沌强度增 加，但对于较大的角动量，混沌变弱。对于带电粒子随动力学参数的变化而引 起的动力学跃迁所得到的数值结果也给出了定性解释。 关键词：天体力学：辛算法；混沌；行星磁气圈 a b s t r a c t t ot r u l yd e s c r i b et h ec h a o t i cm o t i o no f t h eh a m i l t o n i a ns y s t e m ，w en e e dr e l i a 【b l ea n u m e r i c a lm e t h o da sw e l la sc h a o si n d i c a t o r s i nt h i s p a p e r , t h ed e v e l o p m e n ta n d u t i l i z a t i o no ft h en u m e r i c a lm e t h o da r ec o n s i d e r e d p a r t i c u l a r l y a c c o r d i n gt oa l li d e ao ff o u r t h - o r d e rf o r c eg r a d i e n ts y m p l e c t i ca l g o r i t h m sp r o p o s e d b yc h i na n dc o w o r k e r s ，w ec o n s t r u c tt w ot y p e so fe x p l i c i tf o u r t h - o r d e rf o r c eg r a d i e n t s y m p l e c t i ca l g o r i t h m s t h e ya r ee x t e n d e dt os o l v ea l ln b o d yh a m i l t o n i a np r o b l e m w i t hap e r t u r b a t i o nd e c o m p o s i t i o nf o r mo fw i s d o ma n dh o l m a ni nt h es o l a rs v s t e m d y n a m i c s a l s o ，o n eo ft h ed e v e l o p e df o r c eg r a d i e n ts y m p l e c t i ci n t e 霉a t o r si sc h o s e n t os t u d yt h ec h a o t i cm o t i o no fa c h a r g e dp a r t i c l ei nap l a n e t a r ym a g n e t o s p h e r e s o m e d e t a i l sa r ea sf o l l o w s f i r s t ，b ya d d i n gf o r c eg r a d i e n to p e r a t o r st os y m m e t r i cc o m p o s i t i o n s 、v eb u i l da s e t o fe x p l i c i tf o u r t h o r d e rf o r c eg r a d i e n ts y m p l e c t i ca l g o r i t h m s ，i n c l u d i n gt h o s eo f c h i n a n dc o w o r k e r s ，f o ra s e p a r a b l eh a m i l t o n i a ns y s t e mw i t hq u a d r a t i ck i n e t i ce n e r g y 丁 a n dp o t e n t i a l e n e r g yv t h e ya r e h a m i l t o n i a ns y s t e mt h a tc a nb es p l i t e x t e n d e dt os o l v eag r a v i t a t i o n a l n b o d y i n t oak e p l e r i a np a r t h oa n df lp e r t u r b a t i o n p a r tei nj a c o b ic o o r d i n a t e s i ti sf o u n dt h a tt h ea c c u r a c yo f e a c hg r a d i e n ts c h e m e l s g r e a t l ys u p e n o rt o t h a to ft h es t a n d a r df o u r t h o r d e r f o r e s t r u t h s y m p l e c t i c i n t e g r a t o ri nt + v - t y p eh a m i l t o n i a nd e c o m p o s i t i o n , b u t t h e y a r eb o t h a l m o s t e q u i v a l e n tmt h em e a nl o n g i t u d ea n dt h er e l a t i v e p o s i t i o n f o r h o + 日1 由p e d e c o m p o s i t i o n a tt h es a m et i m e ，t h e r ea r en ot y p i c a ld i f f e r e n c e s b e t 、v e e nm e n u m e r i c a lp e r f o r m a n c e so ft h e s eg r a d i e n ta l g o r i t h m s ，e i t h e ri nt h e s p l i t t i n go ft + v o ri nt h es p l i r i n go f h o + 日1 i np a r t i c u l a r , c o m p a r e dw i t ht h ef o r m e rd e c o m p o s i t i o n t h el a t t e rc a nd r a m a t i c a l l yi m p r o v et h en u m e r i c a la c c u r a c y b e c a u s et h i s e x t e n s i o n p r o v i d e saf a s ta n dh i g h p r e c i s i o nm e t h o dt os i m u l a t ev a r i o u so r b i t a lm o t i o n so f n 。b o d yp r o b l e m s ，i ti sw o r t hr e c o m m e n d i n gf o rp r a c t i c a lc o m p u t a t i o n s s e c o n d ，w en u m e r i c a l l yi n v e s t i g a t et h ec h a o t i cm o t i o no fac h a r g e dp a n i c l ei na p i a n e 住呵m a g n e t o s p h e r eb ym e a n so fan e w p r o p o s e df o u r t h o r d e rf o r c eg r a d i e n t s y m p l e c t i ca l g o r i t h mw i t ht h e h 0 + 日i 屯p ep e r t u r b a t i o nd e c o m p o s i t i o n t h e i i a b s t r a c t m e n t i o n e dp h y s i c a lm o d e lc a nb es i m p l i f i e da sa p e r t u r b a t i o nt w o b o d yh a m i l t o n i a n p r o b l e m o nt h ee q u a t o r i a lp l a n e ，s e v e r a lk i n d so ft h ep h a s ep o r t r a i t so ft h ec h a r g e d p a r t i c l es h o u l db ed e t e r m i n e db yt w od y n a m i c a lp a r a m e t e r so ft h er a t i oc h a r g et o m a s sa n dt h ez - c o m p o n e n to ft h ea n g u l a rm o m e n t u m c e r t a i n l y , b e s i d e st h et w o p a r a m e t e r s ，t h ee n e r g ya l s oy i e l d sa ni n f l u e n c eo nt h em o t i o no ft h ep a r t i c l eo u tt h e e q u a t o r i a lp l a n e i ti sf o u n dt h a ti n c r e a s i n gt h ee n e r g yo rt h ea b s o l u t ev a l u eo ft h e r a t i od o e sa l w a y sc a u s et h ee x t e n to fc h a o s h o w e v e r ，c h a o sb e c o m e sw e a k e rf o ra l a r g e rm a g n i t u d eo ft h ea n g u l a rm o m e n t u m s o m eq u a l i t a t i v ei n t e r p r e t a t i o n st ot h e d y n a m i c a lt r a n s i t i o no ft h ec h a r g e dp a r t i c l ew i t ht h ev a r i a t i o no ft h ed y n a m i c a l p a r a m e t e r sa r eg i v e n ，t o o k e yw o r d s ：c e l e s t i a lm e c h a n i c s ：s y m p l e c t i ci n t e g r a t o r ；c h a o s ；p l a n e t a r ym a g n e t o s p h e r e s i l l 目录 目录 第1 章绪论1 1 1 研究历史和现状1 1 2 本文的主要内容和创新点4 第2 章辛算法简介6 2 1 引言6 2 2 非力梯度辛算法6 2 3w i s d o m 。h o l m a n 分解8 2 4 力梯度辛算法1 1 第3 章几个四阶力梯度辛算法1 4 3 1 引言1 4 3 2 四阶力梯度辛算法求解n 体问题1 4 3 2 1 在丁+ y 分解形式下构造四阶力梯度辛算法1 5 3 2 2 推广算法求解摄动k e p l e r 问题18 3 2 3 扩展算法求解j a c o b i 坐标下的n 体问题2 2 3 3 数值模拟2 4 3 3 1 摄动二体问题2 4 3 3 2 三体问题 2 7 3 4 结论和讨论。3 2 第4 章行星磁气圈内尘埃粒子的混沌运动3 3 4 1 引言3 3 4 2 物理模型和数值方法3 3 4 2 1 物理模型3 4 4 2 2s a 4 型辛算法的能量相对误差3 5 4 3 数值模拟结果和诠释3 7 4 3 1 数值模拟3 7 4 - 3 2 定性解释。4 0 i v 目录 4 4 本章小节4 2 第5 章结论与展望4 3 5 1 结论4 3 5 2 展望4 4 致谢4 5 参考文献4 6 攻读学位期间的研究成果4 8 v 第1 章绪论 第1 章绪论 对于任何一个已知其哈密顿方程的物理模型而言，若给出初始条件可以通 过数值方法来求解出系统的各运动量，例如运用传统的数值方法一四阶龙格库 塔方法【l 】。但由于数值模拟存在计算误差，导致运动积分在整个数解过程中得不 到很好的保持。为了能更好的揭示系统的动力学性质，学者对数值计算方法进 行了不断的改进和创新。本文主要工作就是进一步改进数值方法，并将其运用 于研究哈密顿系统的动力学性质。 1 1 研究历史和现状 随着计算机的问世，数值方法得到极大地发展。学者在研究行星动力学引 力n 体哈密顿的数值解问题时，主要是从两个领域着手：一个是用流形改正方 法来处理( n a c o z y 2 l ；f u k u s h i m a 3 】) ，另一个是用辛算法来处( w i s d o m t 4 5 】； w i s d o m & h o l m a n 6 】：w u 【7 , 8 ；z h u i g l ；m c n e i l & n e l s o n t l o 】) 。 为了有效地抑制各种快速增长的数值误差，流形改正方法的机理是在每一 步积分之后，使其得到的数值解乘上一个标度因子来迫使该数值路径回到初始 的超曲面上。对于传统数值方法，例如龙格库塔算法，由于人工耗散或是局部 截断误差的积累所导致的缺陷，当对该算法执行流形改正方法后这些缺陷就不 可能再出现。值得指出的是，用了流形改正的低阶传统数值方法可以被视为一 种快速和高精确度的方案来模拟各种各样的轨道运动( w u 【l l 1 2 1 ；m a 【1 3 1 5 1 ；z h o n g & w u 【1 6 j ) 。 同样，在有效地抑制各种快速增长的数值误差上，辛算法也能保持能量积 分，但不同于前者。它则是保持相空间的辛结构，即所谓的保持结构的几何积 分算法。由于具有这些优势，辛算法特别适合于定性研究天文动力学中的哈密 顿系统长期演化问题。但是辛算法的右函数较复杂，在数值计算时要花费大量 的时间。这将导致其计算效率降低。所以在构造辛算法时不能阶数过高，即精 度达不到很高。这是辛算法存在的缺点。 辛算法的先驱工作可以在八十年代初的文献( r u t h 【1 7 l ；w i s d o m 4 , 5 1 ；f e n g 1 8 1 ) 中 找到。求解町分离哈密顿系统的显式辛算法主要涉及映射方法。作为典型的例 第l 章绪论 子，例女 r u t h t l 7 】将哈密顿分解为动能丁和势能v 两个可积部分的二阶和三阶辛算 法和w i s d o m 【4 5 】在研究与木星3 ：1 共振的小行星的运动中用到的平均原则。然后 有c a n d y & r o z i i l u s 【1 9 】和f o r e s t & r u t h 2 0 】的四阶辛算法以及y o s h i d a l 2 1 1 的高阶辛 算法。应该强调指出的是w i s d o m h o l m a n 6 】著名的交替前进辛算法。它是通过在 j a c o b 坐标下将引力n 体哈密顿分解为k e p l e r 部分日。和摄动部分h 来克服映射 方法的显著局限性。此外，这些丁+ y 分解形式的算法可以直接运用于日。+ 日。分 解形式。众所周知，同一个辛算法在这两种分解形式下其数值表现是不同的： 在精确度上后者大大优于前者。沿着这个思想方向发展出一些先进的方法并在 太阳系问题中得到广泛运用。例女1 w i s d o m t 2 2 j 提出的辛校正方法，c h a m b e r s & m u r i s o n l 2 3 】提出的膺高阶辛算法。这些方法共同的特征是它们都源于非力梯度的 标准辛算法。而非力梯度的标准辛算法除了一阶、二阶算法外都用到一些负积 分步长。 另一方面，三阶力梯度显式辛算法是一个独特的辛算法。它在r u t h r 7 】最早 给出非力梯度的标准三阶显式辛算法时同时出现。这种力梯度算法只含有正积 分步长，但是需要通过力来计算力梯度。它在丁+ y 分解形式中得到实现。然而， 很少有研究人员给予承认或是重视它的应用和发展。这是由于在时间可逆的经 典力学问题中没有必要使用全部的正积分步长或无法方便的计算力梯度。相反， 在求解时间不可逆的方程时对正积分步长的需求至关重要，例如求解含时薛定 谔方程。由于这个原因，c h i n 及其合作者( c h i n t 2 4 - 2 7 ；c h i n & k i d w e l l 2 s 】；c h i n & c h e n 2 9 , 3 0 1 ) 发展了这种力梯度辛算法并把它扩展到四阶。这种算法已经被证实在 精度上比非力梯度的标准同阶f o r e s t r u t h 方法好1 0 至l j 8 0 倍( c h i n 2 4 1 ) ，并且可以非 常有效地解决圆形限制性三体问题。例如求解显时间依赖的力场引力少体问题 ( c h i n & c h e n t 3 0 】) ，以及求解经典力学和量子力学问题( c h i n & c h e n 2 9 1 ) 。本文的 主要工作将在c h i n 的思想基础上继续发展四阶力梯度辛算法，并结合 w i s d o m h o l m a n 的交替前进想法进行推广。 随着数值方法的提出和改进，使得数值研究混沌运动成为可能。混沌运动 是动力天文学中一个重要的研究课题。混沌意味着蝴蝶效应，即最终的动力学 状态是指数敏感于初始条件的微小变动。混沌运动的思想是1 9 6 3 年美国气象学 家l o r e n z 3 1 】对大气对流模型研究时提出。2 0 世纪五六十年代，k o l m o g o r o v 3 2 等 人的研究工作才揭示了保守系统也可能出现混沌，并指出天体运动是确定性和 随机性的对立统一。在天体力学和天文动力学中已经有很多关于混沌的研究， 2 第1 章绪论 这主要是因为混沌在决定动力学结构和太阳系演化问题中起重要作用i j 引。在 w i s d o m 【4 , 5 , 3 4 j 的一系列文章中，用混沌成功的解释了3 ：1 k i r k w o o d 缺口的起源。 s u s s m a n & w i s d o m 3 s 】通过对外行星的8 4 5 百万年的数值积分发现了冥王星的混 沌运动。除了那些复杂的系统外，通常还考虑用来近似代表实际物体运动的两 种简化的动力学模型。一个是平面圆形限制性三体问题，它的相空间结构具有 出奇的复杂度【3 6 】。另一个是处理摄动k e p l e r 问题【3 7 3 8 1 ，由于摄动的存在，在大多 数情况下系统从可积系统变为近可积系统，甚至是高度不可积系统。 作为一个重点，影响混沌鉴定的两个要点是值得一提的。第一，要选择一 个可靠的数值方法。上述两种领域的数值方法都是很好的选择。第二，采用适 当的方法来识别混沌。目前已有一些混沌指标。应该指出的是每个指标都有它 的优点和缺点。例如l y a p u n o v 指数作为一种常见的混沌识别指标，在使用时不 需要考虑系统的相空间维数，但也存在一些优缺点【3 9 喇】。它是随时间变化的两 轨道之间距离与时间的比值，其中最经典的计算方法有两粒子法和变分法。但 这两种经典l y a p u n o v 指数方法需要的重整化次数很多，造成了计算速度慢而且 可能带来更多的计算误差。f r o e s c h l e l 4 3 1 等人提出了三种不需要重整化的快速 l y a p u n o v 指标。其与经典l y a p u n o v 指数相比，它们在速度和性能上都有很大的 优势而且能够区分弱混沌，但对于每次不同的积分初值此类方法都要求解n 次 变分方程。于是f r o e s c h l e 和l e g a t 删对快速l y a p u n o v 指标做出了进一步改进， 使这种方法不仅可以识别混沌和有序轨道，还可以分辨非共振和共振轨道。w ux 4 0 l 等人提出了两附近轨道l y a p u n o v 指数，这种方法不需要求解变分方程，虽然 需要重整化，但是重整化的次数不多，所以计算速度也很快。当相空间的维数 减去独立积分个数不超过三时，p o i n c a r e 截面方法【4 5 】可以为k a m 环或是随机层 的描述提供一个直观清晰的相图。其做法如下：首先在相空间选择一个固定的 截面，如z = 0 截面，让系统的轨迹线沿同一方向多次穿越截面，并记下穿越时 留下的点。分析p o i n c a r e 截面上这些点的分布情况就能了解系统的动力学性质。 当截面上只得到一个点时，说明系统是周期运动：当截面上得到的点构成一条 或是多条封闭曲线时，表明系统做近周期运动；当截面上得到的点是杂乱无章 的排列时，说明系统的运动是混沌的。此外，识别混沌的方法还有很多，例如： 较小排列指标4 6 , 4 7 】、频谱分析【4 8 。5 0 1 、o 1 指标 5 1 , 5 2 1 等等，请见伍【4 5 1 。 第1 章绪论 1 2 本文的主要内容和创新点 本文介绍了辛算法的基本概念和构造格式，并将构造格式运用在c h i n t 2 4 】提 出的四阶力梯度辛算法上，创造了一些新系数的丁+ y 分解形式的四阶力梯度辛 算法。然后把它们推广到由w i s d o m h o l m a n 提出的在j a c o b i 坐标下的日。+ 日。分 解形式来求解引力n 体哈密顿问题。并选用其中一个日。+ 日分解形式的辛算法 来研究行星磁气圈内带电粒子的动力学性质。下面对本文的结构安排做详细说 明。 第二章简要介绍了辛算法的有关概念。较详细的介绍了各种显式辛算法的 具体表达形式及其对应的差分格式。 在第三章中对于分解为二次动能r 和势能v 的可分解的哈密顿系统，通过在 对称组合中加入力梯度算子，我们建立了一组显式四阶力梯度辛算法，其中包 括c h i n 及其合作人的力梯度辛算法。然后这些算法被推广来求解在j a c o b i 坐标 下分解成k e p l e r 部分h 。和摄动部分日。的引力n 体哈密顿系统。通过使用它们 来求解摄动k e p l e r 问题和太阳系n 体问题的数值实验，可以发现在丁+ y 型哈 密顿分解下，每个力梯度算法的精确度明显优越于f o r e s t r u t h 的非力梯度标准 四阶辛方法；但是对于矾+ 日，型分解，在平经度和相对位置精度上它们两个几 乎是等效的。同时，无论是在丁+ v 型分解还是在日。+ 日。型分解中，这些力梯 度算法在数值性能上没有典型的差异。尤其是，后种分解与前种分解相比能大 大地提高数值精确度。 在第四章中根据几种赤道平面的相位图取决于电荷质量比和z 方向角动量 这两个动力学参数的情况，我们数值研究了行星磁气圈内的带电粒子的混沌运 动。我们的目标是应用上一章提出的s a 4 型辛算法通过研究电荷质量比，角动 量和能量这三个参数的一些可能值来数值描绘赤道平面外的轨道动力学。本系 统的混沌主要依赖于这三个参数中的任何一个。又由于本系统的对称性，系统 自由度从三维减少到两维。在这种情况下，我们选用p o i n c a r e 截面方法来描绘 本系统的混沌运动。数值实验发现当增大能量或是比率的绝对值时，总是导致 一定程度的混沌。然而，对于大量级的角动量，混沌是较弱的。并对得到的结 果给出了定性解释。 第五章是对本文工作的总结和对以后工作的展望。对于在对称组合中加入 力梯度算子的显式四阶力梯度辛算法，我们已经做了比较全面的推导，并将它 们推广到哈密顿摄动分解体系。然后选用它们来研究行星磁气圈内的带电粒子 4 第1 章绪论 混沌运动。未来，我们将更加关注辛算法的发展，提高对辛算法、混沌识别方 法的认识，努力做出创新。 本文的创新点： ( 1 ) 根据辛算法的构造格式，推导完善了在对称组合中加入力梯度算子的 r + y 分解形式的四阶力梯度辛算法。并将其推广到求解可以在j a c o b i 坐标下分 解成k e p l e r 部分h 。和摄动部分h 。的引力n 体哈密顿系统。 ( 2 ) 将占一a 4 型辛算法与p o i n c a r e 截面结合起来，对赤道平面外的行星磁气 圈内带电尘埃粒子模型进行数值求解。并多次改变各物理参数的取值，根据 p o i n c a r e 截面上点的分布情况来判断各物理参数对该系统混沌性质的影响。 第2 章辛算法简介 第2 章辛算法简介 2 1 引言 对哈密顿系统混沌运动的描述，有赖于可靠的数值方法和识别混沌的方法。 正如绪论部分指出，传统数值方法由于人工耗散等因素导致运动积分在整个数 解过程中得不到很好的保持，而辛算法特别适合于定性研究天文动力学中的哈 密顿系统长期演化问题。 辛算法是由f e n g t l 8 1 和r u t h 1 7 1 为了克服传统算法不足而各自独立建立的算 法。这种算法能保持住哈密顿系统的相空间辛结构，即所谓的保持结构的几何 积分算法，简称辛算法。前者主要运用正则变换的隐式生成函数从事隐式辛算 法的研究工作。隐式辛算法适用于任何形式的哈密顿函数，但由于在计算中用 到了迭代方法，而使得它的计算效率降低。后者率先对哈密顿函数能分解为动 能r 和势能矿的情况提出了显式辛算法。由于显式辛算法的计算效率比隐式辛算 法的高，当哈密顿系统能可积分解时，显式辛算法是更好的选择。本章将主要 介绍显式辛算法的具体表达形式及其分类。 2 2 非力梯度辛算法 r u t h t l 7 】对哈密顿函数能分解为动能丁和势能矿两个可积部分的情况提出的 显式辛算法形式如下： e c e = e 彤 、 e r v e 仃= e r w ( 2 2 ) 这是最简单的丁+ y 分解形式的一阶显式辛算法，其中于和矿是l i e 算子。式( 2 1 ) 与式( 2 2 ) 的区别在于于算子、矿算子的排列顺序不同。这将导致其差分格式不同， 即计算次序不同，详见式( 2 5 ) 和式( 2 6 ) 。 对于具有如下分解形式的哈密顿 h = r ( 矽) + y ( 虿)( 2 3 ) 假设丁只是有关动量卢的正定二次函数，例如丁( 多) = 多2 2 ；v 代表一个位置函 6 第2 章辛算法简介 数牙。我们也定义l i e 算子于和矿分别与r 和v 相关，其中于= d r = ，n ， 矿= d v = ，n ( 符号 ，) 是泊松括号) ，相当于， 于= 喜p 云，矿= 喜瑶 亿4 ， 其中，= - 0 z 为，q 。和p 。分别表示第，个广义坐标和第f 个广义动量，厂是力 的第f 个组成部分。由正则方程我们可以很容易写出式( 2 3 ) 在( 2 1 ) 和( 2 2 ) 形式下 的差分格式 霞= 磊一。+ 矿( 玩一。) 玩= 玩，+ 或( 2 5 ) 和 玩= 玩一l + 识一l 庶= 觅一l + 矿( 玩) ( 2 6 ) 其中f 为积分步长，s ( - - 1 , 2 ，) 是积分步骤的标号。而传统的e u l e r 数值方法其差 分格式如下 夙= p 川+ 矿( 玩一。) 玩= 玩一l + 破一l ( 2 7 ) 通过比较得出：前种数值方法的两个积分是相关联的，即第一个积分的结果是 第二个积分的初始值；而后种传统e u l e r 数值方法的两个积分是相互独立的。并 且可以证明，前者的外积满足氟a 姬= 纸一。 氟一。这说明式( 2 1 ) 和式( 2 2 ) 保持了相空间的辛结构。 丁+ 矿分解形式的二阶显式辛算法 p 也r ；e j l 矛= p r w ( 2 8 ) 和 p 三f l ，p 开p j lf l ，= p r w ( 2 9 ) 是在式( 2 1 ) 、( 2 2 ) 的基础上运用映射方法得到。式( 2 8 ) 、( 2 9 ) 差分格式分别如下 一。一 l q 2 q s - l + j 矾一l 7 第2 章辛算法简介 p ，= 虎一。+ f f ( q + )( 2 1 0 ) 。1 吼2 口+ i 识 和 一：ps-i+pp 委矿( 磊一。) 2 i 秒( 吼一1 ) 玩= 玩一l + r p 。( 2 1 1 ) 庶：p 一+ 丢矿( 玩) 只2 + j 可( 吼) 高阶辛算法阵1 y o s h i d a 2 1 1 仍在此基础上连续使用对称方法得到。由于它们都是式 ( 2 1 ) 、( 2 2 ) 的映射，所以它们都能保持相空间的辛结构。而这些非力梯度辛算 法除了一阶、二阶算法外都用到一些负积分步长。这在求解时间不可逆问题时 会遇到应用的困难。 2 3w i s d o m h o l m a n 分解 w i s d o m h o l m a n 6 】的著名辛算法是通过在j a c o b i 坐标下将引力n 体哈密顿 分解为k e p l e r 部分日。和摄动部分日。来克服映射方法的显著局限性。其分解形 式如下 h = h o ( 虿，多) + s h l ( 虿) ，占 1( 2 1 2 ) 这里占是作为一个让日远小于日。的小的摄动参数。众所周知，哈密顿函数的不 同分解方式会影响辛算法的精度。例如，上节提到的r u t h 1 7 1 的丁+ y 分解形式的 n 阶非力梯度辛算法的解误差为o ( r 胂1 ) ；而当哈密顿采用w i s d o m h o l m a n 分解 形式( 2 1 2 ) 时，其辛算法的解误差为o ( s r 1 ) 。通过比较可以发现在同一阶数n 下，后者的精度和效率都优越于前者。所以w i s d o m h o l m a n 的这种对太阳系n 体问题所给出的哈密顿函数摄动分解情形是值得引入辛算法来提高其精度的。 w i s d o m h o l m a n 的哈密顿分解方式详细介绍如下： 在质心坐标系下，太阳系n 体问题的哈密顿函数为 h = 芝虹2 一宝窆竺生生(213)m 儡 f- i = l 等lq f ， r 7 第2 章辛算法简介 其中q 口爿玩一牙l 是天体f 和天体之间的距离；江0 表示太阳，m ，是天体f 的 质量；g 是引力常数。鉴于这种情况，w i s d o m h o l m a l l 【6 j 给出的方法成功的实施 了所需的分解形式。通过利用雅可比坐标有 牙 玩一杰矽“硼) ( 2 。4 ) 赤善聊厄 o 0 ) 其中 m = 肌 ( 2 1 5 ) 一 j 、 ，z ；= - 1 m “, m 0 忙 i 2 善3 ( p 。砉+ 晶旦a p , ) ，直= ，日。 = 善3 ，瓦3 2 1 ) 其中g ，= - a r o a q 。，z = - a h , a q ，。如果上述的于和矿分别被n o 和丘取代， r + y 分解形式的力梯度辛算法可以很容易的重新写出它的表达形式。例如，下 面我们给出算法彳l 和b i 的修改版本 沪卸p 吾p 昙硝p 玑三她+ 面t 3 竹砒巾p ；p 言豳p 吾成如) 和 p b 1 ) p 吾以+ 未占2 【扁，成，加p 詈戌p 言矾p ；吼p 言矾p ；成p a + 面t 3 占2 【疗，成，疗- 】 ( 3 2 3 ) 使用g 一型算法的关键点是如何处理 疗，疗。，疗。】算子。事实上，这个算子的计算 和【矿，于，矿】算子的计算方法是完全相同的，也就是说， q ，风，q 】= ， q ，h o ，q ) = ， q ，t + v o ，日。) ) “q 刀，q ) = 2 善3 瑶云= 善3v ，旧2 吾( 3 “) 3 2 2 2 离散差分格式 根据哈密顿函数的不同分解方式，我们从两种观点设计了这些力梯度辛算 法的差分格式。 情况1 ：丁+ y 分解。对于分解为动能丁和势能y = t o ( q ) + z h 。( 虿) 的摄 动k e p l e r l h 题，我们有z = 一a v 却，。然后，我们仍然用算法4 l 和b 1 为例来说 明差分格式的机制。等式( 3 1 2 ) 的差分格式如下 ( 彳1 ) q5 q ，一1 + 7 p 。一1 n 1 9 第3 章几个四阶力梯度辛算法 一2os-ipp1 詈矿( 虿) 2 i 矿【gj 一一r 一 口2g + = p j p 一= 卢+ 吾t i n ) + 石i f 2v 厕) 1 2 】 ( 3 2 5 ) 一一f 一 g 2 g+ = p j p ，= = p 。- + 言矿( ；f ) ，2 + i 矿l gj g s2g + i 几 其中s ( = 1 , 2 ，) 是积分步骤的标号。同时，等式( 3 1 9 ) 的差分格式如下 ( b 1 ) p 一= p s - i + 如+ 毛v i 加汗】 一一f 一 g 2 q ，一i + = p j p + + = = p + 言矿( 虿+ ) p 2 p + = 可l gj q j = 一+ j (326)q p2 6 ) 2 + ： 【 p = = p + 言r j j ：( 虿) 2 + = 旧j 一一f 一 qs 2 q 七ip j 多，= 矽+ 争地) + 石t 2 v7 ( o 。) 门 情7 况u 2 h 。+ h 。分解。当哈密顿有一个像等式( 3 2 0 ) 那样的可分离的形 式时，为主的k e p l e r 部分h 。从初始状态( 磊，风) 经过时间丁后有分析解。这个分 析解用高斯函数表示为：虿= k e p l e r l ( f 1 0 ，风，f ) 和卢= k e p l e r 2 ( # o ，死，r ) ，其具体 解法见2 3 节的式( 2 2 0 卜弋2 2 2 ) 。考虑到夕= 二御。衡，7 y 程( 3 2 2 ) 和( 3 2 3 ) 符合以 第3 章几个四阶力梯度辛算法 下差分格式： 和 多”= 矿+ 百t e 晰- ”) + 百t 2 c v i 加) | 2 】 ( 3 2 7 ) 多= 广+ 詈巧( n 虿，= k e p l e ，1 ( 牙”，p + + + ，吾) = ，k e p l e ，2 ( n p + 4 - - i - ，f ( 占一b 1 ) p 。咆一- i - t f ( o 一) + 害v i 厕汗】 2 l 一阼 一艮 一旧 一 一 刊 ， 飞 。 。 d 卜尹 嘶 蜥 一坤一一 = = 一p ， = 。= 一p g p 艮 切 。 脚 脚 = = 一g j p 、， 苁一y o 多， 萨 弘 玩 巾 孙 尼 盯 鳓 础 啦 却 矿 矿 第3 章几个四阶力梯度辛算法 p = p + 詈巧( 虿) 多“= k e p l e ，2 ( 虿，多。，要) 卢一= 多”+ 昙亏( 虿”) 牙，：k e p l e ，1 ( 矿，p 嘲，三) p 。= k e p l p ，2 ( 矿，多。，要) 驴多+ 詈魄) + 等v i 加川2 】 h ：窆世2 一窆守竺竺生(329)m 盏 f i = 1 一l - - + , l q 盯 r 7 弘摩等 删 1 一般而育，对于等式( 3 2 9 ) 这种哈密顿分解形式若直接在质心坐标系下积分是不理想的。这足冈为太阳过 于接近质心。相反，数值积分通常在f 心系内实施，这时太阳的演化就町以小考虑。在每一个积分步后， 太阳的演化足通过质心关系和总动量守恒来考虑的。 第3 章几个四阶力梯度辛算法 其中 m ，= 聊， m ：= 心- i m 川i m o 忙 i 。) n ( 3 3 1 ) ( 3 3 2 ) 他1 i j 酮嗡苗坝里耕与力 h = h 0 + 凰 ( 3 3 3 ) 凰= 善c 等一半， 3 4 ， 即艺( g m r o m , 一警，一喜善警 3 5 ， 说明：q i 虿川和当质心坐标虿，在等式( 3 3 0 ) 中被表示成j a c o b i 坐标虿：后，它就没 有再出现在等式( 3 3 5 ) q h 。h 。是k e p l e r 部分，h 。表示行星间的相互作用部分。 值得指出的是前者远远大于后者。 通过对方程( 3 3 3 ) 和( 3 2 0 ) 1 拘比较，可以很容易的理解在j a c o b i 坐标下力梯度 辛算法也适用于n 体系统。下面是对这个问题的简要介绍。 对于在j a c o b i 坐标下的t + v 分解，动能t 和势能y 描述为 丁= 莩等 ( 3 3 6 ) 丁= 等 ( 3 y ：一芝鱼型生一芝艺g m i m j ( 3 3 7 ) - 1 = 1 q j o智等lq “ 、。 然后，我们得到 于= 莩蔷3 鲁p v 可0 ，矿= 善n - i 缶3 厶两 b 3 8 ， 眈杪m t 艺j = lj 童, k = 鲁嚣毒= 荟n - i 蔷3v 口c 莩警，蠢c 3 舯， 在上面的等式中：= 一a 吖却：；指标f ( 或者，) 代表行星f ；指标_ ，( 或者后) 表示第个部分；符号v 。是梯度；第，个坐标由行星，组成。在这种情况下，无 第3 章几个四阶力梯度辛算法 论是像等式( 3 1 2 ) 的a 型算法或是像等式( 3 1 9 ) 的b 型算法都可以像运用于摄动二 体问题那样容易的实施。 另一方面，对于等式( 3 3 3 ) 习g 样分解的哈密顿，实施像等式( 3 2 2 ) n ( 3 2 3 ) 另g 样的型算法是典型的类似于摄动二体问题