积分器和程序设计.doc

紫金山天文台台刊稿件 PMOE 2003行星历表框架 （II）积分器和程序设计 李广宇 田兰兰 （ 中国科学院紫金山天文台，210008） 摘要: 积分器是实现历表数学模型的必要工具. 程序设计则是实现这两者，进行历表研究和完成实用历表的最终手段. 本文从历表研究需要和理论与实践结合的角度，统一整理了有关的知识和方法，详细介绍了张家祥积分器的理论基础，推导了实用公式，给出了有关数据结构、算法和程序设计的原则与方法，并对积分器的精度进行了分析. 关键词：积分器，数据结构，程序设计 数值算法 我们在文[1]中叙述了PMOE 2003精密行星历表框架的数学模型，对框架所用的积分器和程序结构及设计问题，只作了概要介绍. 作为历表框架的重要组成部分，积分器是实现历表数学模型的必要工具. 程序设计则是实现这两者，进行历表研究和完成实用历表的最终手段. 积分器和程序设计的质量直接影响历表的实现和质量. 有关这两部分内容的叙述散见于天体力学[2-5] 、计算方法[6,7]和程序设计[8-10]的文献中. 张家祥先生长期从事太阳系动力学模型数值积分的研究，积累了处理各种问题的丰富经验. 他改进Cowell方法建立的积分器具有快速高效和大幅减少舍入误差的显著优点，我们已在文[1]的引言中概要介绍了他的有关研究工作. 张先生关于积分器的工作在文[11,12]中已有介绍，但对学习者尚觉简略，并且没有涉及程序设计问题. 近年来作者与倪维斗等合作研究精密行星历表框架，移植了张先生的积分器程序库[13]，并在面向对象和结构化程序设计的原则下进行了改进，完成了Pascal[14]和C++[15]两种语言的版本. 在此过程中对数值积分和程序设计的有关问题进行了整理，其间多次与张先生讨论，并参阅了先生的有关笔记. 本文试图从历表研究的需要和理论与实践结合的角度对有关的知识和方法进行统一的整理和陈述，以满足学习者的需要. 设 q=(q1, q2,…, qN)τ为天体坐标，天体的运动方程一般具有如下形式[2]： （1） 本文讨论用Cowell 数值积分方法[2.3]求解这类方程的问题. 数值积分的对象是单个坐标分量，所需要讨论的，只是积分方程组(1)中的一个方程的方法.因此我们略去坐标分量和右函数的下标i，将方程写为的形式以求简明. 本文第一部分用符号算法[2,3,6]导出构造张家祥积分器必需的公式，第二部分给出用于积分器程序设计的数据结构和算法[8-10]，最后对积分器的精度加以讨论. 有限差算子[3,6] 数值算法依赖于称为步长的实数h，命函数 （2） 并假定f(t)具有一切必需的分析性质. 定义如下4个作用于f(t) 的有限差算子： 移位算子 E f(t)=f(t+h) 中差算子 δf(t)=f(t+h/2)-f(t-h/2) 平均算子 μf(t)=[f(t+h/2)+f(t-h/2)]/2 不变算子 I f(t)= f(t) 加上熟知的微分算子 , 一共是5个算子. 由定义可以得到联系微分算子和移位算子的关系： E=ehD （3） 由此推论：对任一实数m有 按定义， 故有 ，两式相乘得 ，解出μ即得联系中差算子和平均算子的关系： （4） 又 ，平方之，并以(4)式代入，即得联系中差算子和移位算子的关系： （5） 由(4)式可展开平均算子μ和μ－1为中差算子δ的幂级数： （6） （7） 差分 给定实数t0，以h为步长，取一列等间距步点： 记函数 f(t) 在这些步点上的对应值为,以中差表示的此函数的差分表为： （8） 记号，余类推. 显而易见，表中各元素间有如下的关系：(i) 任一元素为其左下方与左上方两元素之差； (ii) 任一元素为其上方和右上方两元素之和. 表8还需要向左扩展出和两列. 给定一个常数做为元素的值，按上述关系 (ii) 即可构造出列的其余元素. 同样，再给定一个常数做为元素的值，则可以把上表再向左扩展一列，构造出列的其余元素. 和列的元素通常分别称为函数 f(t) 的一阶与2阶和分，并分别记为和. 为简便起见，用统一的记号表示函数及其差分为： (9) 函数符号的上标表示差分的阶数，下标表示对应步