积分器和程序设计.doc

Cowell 数值积分方法紫金山天文台台刊稿件PMOE 2003行星历表框架（II）积分器和程序设计李广宇 田兰兰 （ 中国科学院紫金山天文台，210008）摘要: 积分器是实现历表数学模型的必要工具. 程序设计则是实现这两者，进行历表研究和完成实用历表的最终手段. 本文从历表研究需要和理论与实践结合的角度，统一整理了有关的知识和方法，详细介绍了张家祥积分器的理论基础，推导了实用公式，给出了有关数据结构、算法和程序设计的原则与方法，并对积分器的精度进行了分析.关键词：积分器，数据结构，程序设计1 数值算法我们在文1中叙述了PMOE 2003精密行星历表框架的数学模型，对框架所用的积分器和程序结构及设计问题，只作了概要介绍. 作为历表框架的重要组成部分，积分器是实现历表数学模型的必要工具. 程序设计则是实现这两者，进行历表研究和完成实用历表的最终手段. 积分器和程序设计的质量直接影响历表的实现和质量. 有关这两部分内容的叙述散见于天体力学2-5 、计算方法6,7和程序设计8-10的文献中. 张家祥先生长期从事太阳系动力学模型数值积分的研究，积累了处理各种问题的丰富经验. 他改进Cowell方法建立的积分器具有快速高效和大幅减少舍入误差的显著优点，我们已在文1的引言中概要介绍了他的有关研究工作. 张先生关于积分器的工作在文11,12中已有介绍，但对学习者尚觉简略，并且没有涉及程序设计问题. 近年来作者与倪维斗等合作研究精密行星历表框架，移植了张先生的积分器程序库13，并在面向对象和结构化程序设计的原则下进行了改进，完成了Pascal14和C+15两种语言的版本. 在此过程中对数值积分和程序设计的有关问题进行了整理，其间多次与张先生讨论，并参阅了先生的有关笔记. 本文试图从历表研究的需要和理论与实践结合的角度对有关的知识和方法进行统一的整理和陈述，以满足学习者的需要. 设 q=(q1, q2, qN)为天体坐标，天体的运动方程一般具有如下形式2： （1）本文讨论用Cowell 数值积分方法2.3求解这类方程的问题. 数值积分的对象是单个坐标分量，所需要讨论的，只是积分方程组(1)中的一个方程的方法.因此我们略去坐标分量和右函数的下标i，将方程写为的形式以求简明. 本文第一部分用符号算法2,3,6导出构造张家祥积分器必需的公式，第二部分给出用于积分器程序设计的数据结构和算法8-10，最后对积分器的精度加以讨论.1.1 有限差算子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）1.2 差分给定实数t0，以h为步长，取一列等间距步点：记函数 f(t) 在这些步点上的对应值为,以中差表示的此函数的差分表为： （8）记号，余类推. 显而易见，表中各元素间有如下的关系：(i) 任一元素为其左下方与左上方两元素之差； (ii) 任一元素为其上方和右上方两元素之和.表8还需要向左扩展出和两列. 给定一个常数做为元素的值，按上述关系 (ii) 即可构造出列的其余元素. 同样，再给定一个常数做为元素的值，则可以把上表再向左扩展一列，构造出列的其余元素. 和列的元素通常分别称为函数 f(t) 的一阶与2阶和分，并分别记为和. 为简便起见，用统一的记号表示函数及其差分为： (9) 函数符号的上标表示差分的阶数，下标表示对应步点的序数. (9)式中的阶数也可以是负整数，对应于函数的和分. 1阶和-2阶差分即是1阶和2阶和分，分别记为和.应该注意，按照定义，步点序数增加的方向取决于步长h的正负，只有当h为正时，序数增加的方向才与自变量t增加的方向相同.由表可见，偶数阶差分都位于整步点上，奇数阶差分都位于半步点上.由于下面的讨论需要在整步点上进行，因此特别定义整步点上的奇数阶差分为相邻两半步点上同阶差分的平均值，并用新记号表示为： (10)易知奇数阶差分和偶数阶差分之间有如下关系： (11) 现在，整步点上的函数值及其各阶差分与和分可以统一写成如下形式： （12） 从左到右的各列依次为二次和分、一次和分、函数值和一阶至五阶差分.在式(11)中取k=0, i=0，即得到联系一阶和分与二阶和分的关系： (13)由表(12)各元素间的关系，可以导出如下根据二次和分在第i-1和i两个步点上的值，外推出第i +1步点上的值的公式： (14)易知，只要给定了与两个元素的值，同样可以构造出一阶与二阶和分的其余元素，这正是积分器实际使用的方法.1.3 用差分表示微分由(3)(5)两式有：，由此 ，求导得 , 积分之，即得联系中差算子和微分算子的关系： (15)以被积函数的幂级数表达式（7）代入，即得到以中差算子表示微分算子hD 的展开式： （16）显然hD为的奇次幂级数. 算子hD多次作用时，要注意次数的奇偶. 如果称只有奇次项的幂级数为奇次幂级数，只有偶次项的幂级数为偶次幂级数， 那么两个幂级数相乘时乘积级数的通项为从每个乘数级数中任取一项相乘所得的积，其指数为相乘各项指数之和. 由于奇数个奇数之和为奇数，偶数个奇数之和为偶数，可以得到如下关于幂级数之积的结论：奇数个奇次幂级数之积为奇次幂级数，偶数个奇次幂级数之积为偶次幂级数. 由于为偶数个奇次幂级数之积，为偶次幂级数，其各项为的偶次幂，均在整数步点上. 但为奇数个奇次幂级数之积，为奇次幂级数，其各项为的奇次幂，均在半步点上. 为了移到整步点上，根据(4)式和(10)式，只需以因子 乘之，使各项的具有2k+1的形式.一般地，有 （17） （18）1.4 微分算子的展开式由(16)式可知,如上两个偶数阶和奇数阶的微分算子可以表示为中差算子的奇次幂级数的乘积,从而分别展开为偶次幂级数和奇次幂级数,并且可以进一步用待定系数法定出幂级数各项的系数.由幂级数乘积的性质可知,的展开式中幂次最低的项为i ,其系数为1. 因此展开式可以分别写成如下形式： （19） （20）式中 其余系数待定.由于，故有：对（17）式两边求导可得： 即，以（19）、（20）两式代入得：， 比较两边2i的系数，即得偶数阶算子和相邻奇数阶算子展开式的系数之间的如下递推关系： (21) 类似地，对（18）式两边求导,注意到 得:，以(19)、(20)两式代入得：，注意到 得：，两边以4乘之，注意到得：，比较2i项的系数即得奇数阶算子和相邻偶数阶算子展开式的系数之间的第二个递推关系： (22)由递推式 (22) ,可以由奇数阶算子展开式的第一个系数（由(20)式, 此系数为1）和低一阶偶数阶算子展开式的系数推出其余的系数. 由递推式 (21) ,可以由奇数阶算子展开式的系数推出高一阶的偶数阶展开式的系数.由数学归纳法可以证明，1阶算子展开式的系数： （23）由(21)式立即得出2阶算子展开式的系数： （24）由此可得1、2两阶微分算子的展开式： （25） （26）由递推公式 (21) 和 (22)，可以进一步得到各高阶算子的展开式： (27) 1.5 积分公式按定义， (28) (29)以（23）式和（24）式代入，用长除法即得： (30) (31)以算子和 作用于（2）式给出的函数f(t)，即得到速度和坐标的积分公式： (32) (33)1.6 用函数值表示差分实际应用中给出的往往不是各阶差分，而是函数f 在各步点上的值，因此需要以这些函数值来表示（32）和（33）式中的各阶差分.关于差分与函数值的关系，对于偶数阶差分，可以用数学归纳法证明： （34）易知，展开式的系数与二项式的展开式的系数相同.对于奇数阶差分，由（10）式有： （35）式中约定，当i2k时，由此可得用函数值表示各阶差分的系数如下：奇数阶差分：偶数阶差分：由差分表（8）和（34）式还可导出准确至8阶差分和准确至10阶差分的外推公式： （36） （37） 1.7 整数步点上坐标、速度的表达式给定了天体在初始时刻t0的坐标q0和速度，就可以用数值积分方法求解方程(1)，得到天体在指定时刻t的坐标和速度. Cowell方法属于多步法，求解工作分为两步：对(2)式所给函数，也即天体一个坐标的数值积分分两步进行：第一步，根据初值准备在至共9个步点上的值 和+4、+5两个步点上的二次和分，称之为开头；第二步，从已构成的函数表与和分表出发，逐步外推，求出每一步点（即+5步点）上的坐标速度，并更新函数表与和分表，直至指定时刻t，称之为外推.对于二体运动，各量可以直接计算得到，开头的任务十分简单. 对于多体运动，没有分析公式可以直接计算出天体在各个步点上的坐标、速度、加速度和对应的函数值，需要采用迭代算法. 方法是先用二体运动公式计算各步点的函数值作为近似，再用数值积分迭代计算直至收敛. 迭代的次序是：先从0步点开始，计算出和分初值，逆向递推（由t向t-h的方向）4步至-4步点，再从0步点开始，计算出和分初值，正向递推（由t向t+h的方向）4步至+4步点，如此反复进行.下面推导积分器实用的整数步点上坐标速度的表达式. 在（33）式中取差分至f（6），可得q0 的表达式如下： （38) 向下或向上平移一个步点，即可得q1和 q1的表达式：（38）式向下平移两个步点，可得q2的表达式：以f5的表达式(36)代入得：同法可得q2的表达式：依次推出q3, q4, q5和q-3, q-4, q-5的表达式. 数值积分开头时，坐标初值是给定的，因此（18）式用于计算和分的初值，其余各式则用于计算坐标. 各式可以统一写成11： (39)式中，,， ， (40)在(32)式中取差分至f（7），可得的表达式如下：, 由(11)式有 ，代入得： （41）逆推积分时步长为-h，对应于各步点的下标亦须改变符号： （42）用与推导坐标表达式相同的方法，可以得到：等. 与坐标表达式的情形相同，(41)、(42)两式分别用于在开头时计算和分初值.与. 其余各式则用于计算速度. 各式可统一表示为如下形式11.需要注意的是，与（42）式的情况相同，在数值积分开头逆推迭代时，各式中的步长h应换为-h: (43)式中， (44) 将矩阵V的第一行第5列的元素-1814400改为1814400，即得到矩阵V .在数值积分外推过程中，每一次得到的和分是，而不是，这时需要用前者表达出后者.由(14)式有可以导出： （45） （46）代入第 (38)、(42)式，又注意到 ，遂得： （47） （48） 此二式用于在输出前校正中心步点上的坐标速度数据.1.8 非整数步点上坐标速度的表达式3设实数n满足，本节将导出对应于时刻t的坐标和速度的表达式. 显然，时刻t位于步点i和i+1之间.由(3)式，对实数i + n有： （49） （50）当k=2 时有： (51) 以的表达式(25)-(29)各式代入化简得： (52)应用差分算子的性质(11)式消去奇数阶差分，并命m=1 - n，即得坐标表达式： (53)对t求导，注意到，即得速度表达式： (54)不难看出，整数步点上的表达式(39)、(40)式是上式当n = 0时的特例.取即可得到在中心点附近的坐标和速度表达式，（55） （56） 2 程序设计2.1 数值积分的数据结构从本节起，讨论用函数(2)在9个步点上的值对天体数据进行数值积分的具体算法和积分器的程序设计问题.数据结构是算法和程序结构的基础8，我们先从数据结构开始讨论.程序中需要先说明以下全局变量：变量名类型意义ninteger天体数t0extended起始时刻hextended步长nninteger步点序数nneinteger步点总数 整数步点上坐标和速度表达式系数分子的矩阵U说明为常量数组：和和 const US:array0 . 5, -4 . 4of extended；速度表达式系数分子的矩阵V和V说明为常量数组： const VS:array-1.5,-4.4of extended；为了计算和分初值时编程的方便（参见后面2.2节(ii)），数组VS中下标为-1的行对应于矩阵V的第一行，下标为0的行对应于矩阵V+的第一行，下标为15的各行对应于两个矩阵的第26行.数组的第一下标对应于坐标速度的步点序数，第二下标对应于函数值的步点序数. 公共分母说明为常量COMMONDEN = 3628800.0.处理一个天体的运动所需的数据结构和方法封装在如下的类(class)中： TBody=class(TObject) BName: string20; Style: string10; m: real; x0, v0, x, v a: TVector; Elmn: array0.7 of extended; f: TAccel; sf: TSum; constructor Init (t:TTable); Procedure CalcuIniss (c: integer); procedure CalcuSs (const k: integer); procedure CalcuXV (const c: integer; const ux, uv: array of real); procedure CalcuXV (const c:integer;const ux,uv:array of real); overload;function CalcuXV (const nf: extended; const velosity: boolean):TVector; overload; procedure ReCalcuXV; procedure SendF (const i:integer); procedure ElemToCor (const t:extended); procedure CorToElem (const t:extended); procedure CorToCor (const dt,CenterGm:extended); end;前9个数据分量为天体的基本数据：BName:string20;天体名Style;天体类型mu: real;天体质心引力常数x0, v0,初始坐标、速度x, v, a: TVector;坐标、速度、加速度向量elmn: array 0.8 of extended ;6个轨道根数、偏近点角、真近点角、向径向量类型Tvector说明为有3个分量的一维数组： Tvector = array1 . 3 of extended;第10个数据分量 f: TAccel为函数表，存放函数（2）在9个步点上的值，其类型说明为： Taccel = array 1 . 3, -4 . 4 of extended;任一时刻，数组中都存放着9个步点的数据，第一下标对应于坐标分量的序数，第二下标对应于步点的相对序数，即当前步点相对于中心步点的序数. 第11个数据分量sf:Tsum为和分表，用于存放二次和分.这是一个数组： Tsum = array 1 . 3，1 . 2, 1. 3 of extended;数组的第一下标仍然对应于坐标分量的序数，第二下标对应于步点的相对序数(以下提到步点序数时如无特别说明，均指相对序数).对每一个坐标分量，存放对连续两个步点的二次和分；第三下标对应于每个和分的三个分量：第三分量存放该和分之值；该和分与一个预定义的因数ZHFACTOR之积的整数部分存放于第一分量，小数部分存放于第二分量 .数值积分每一步外推生成和分的过程中，总要将做为小量的函数值的和数加到原来的和分上去，直接相加会损失小量的有效位数，产生不容忽视的舍入误差，这是数值积分累积误差的一个重要来源.为减小和消除这一部分误差，张家祥提出了以下处理：(i)和分乘以整数因数ZHFACTOR后分离整数部分和小数部分分别存入两个分量，以增加有效数字的长度.(ii)在有和分参加的加减运算中，特别注意避免小量和大量的直接加减.先分离大量的整数部分和小数部分，然后再以其小数部分和小量加减. 我们引入上述数据结构，就是为了保留尽可能多的有效数字，延缓舍入误差的积累和传递.下面用于处理和分与其他量相加的过程AddSum的算法即体现出了这种思想.过程在分离一数的整数和小数部分时调用了系统函数Tranc和Fcedure AddToSum(var sum:array of extended;w:extended);var w1:extended; k:integer;begin k:=low(sum); w:=zh*w; w1:=sumk+1+frac(w); sumk:=sumk+trunc(w)+trunc(w1); sumk+1:=frac(w1); sumk+2:=(sumk+sumk+1)/zh;end;和分sum的三个分量以数组参数传递给过程，数组的最小下标k用函数low（）取得.这样的处理显著减小了误差积累（见第3节）.用于数值积分的方法有：名称类型参数功能SendFprocedurei:integer更新函数值表CalcuIniSsprocedurec:integer计算和分表初值CalcuSsprocedurek:integer更新和分表CalcuXvprocedurec:integer;const ux,uv:array of extended计算整数步点上的坐标速度CalcuXvfunctionnf:extended;const velosity:boolean计算非整数步点上的坐标速度CalcuXvprocedure无计算中心步点上的坐标速度名为CalcuXv的方法有3个，需要传递的参数各不相同，分别用于计算天体在整数步点、非整数步点和中心步点上的状态(即坐标和速度). 这是用了函数重载技术.另有有三个方法用于坐标计算16：名称类型参数功能ElemnToCorproceduret:extended由根数计算状态CorToElemnproceduret:extended由状态计算根数CorToCorproceduredt,CenterGm:extended由时刻t的状态计算时刻t+dt的状态所有天体的对象说明为一个动态数组：x:array of Tbody;在程序启动后根据实际天体的数目n为它动态分配内存.2.2 数值积分的主要类方法以下按功能说明用于数值积分的6个类方法：(i) 更新函数值表方法：procedure SendF(const i:integer);功能：参数i为步点序数.以加速度向量a的值更新函数值表f的第i列，即对应于第i步点的函数值.算法：见(2)式.以程序语言表示为：for j:=1 to 3 do fj,i:=h*h*aj; 向量a存放天体在第i步点的加速度，即函数之值.这是调用过程CalcuAccel完成的.有关此过程的算法，在文1中详细介绍.(ii) 计算和分表初值方法：procedure CalcuIniSs(c:integer);功能：根据函数值表计算和分表初值.参数c的取值为1,c=1时，正推积分，计算；c=-1时，逆推积分，计算.算法：计算用（38）式，系数在常量数组US0,k:计算用（41）和（42）式，系数常量数组在VS0,k和VS-1,k:下标表达式借用了程序语言的表达方式，此表达式的值当时为0，当时为-1.算法用Warnier图9表示如下：(iii) 更新和分表方法：procedure CalcuSs(const k:integer);功能：参数k为步点序数，过程根据第k个步点的函数值更新和分表.算法：根据二次和分外推公式（14）：计算，过程入口时sfj,1和sfj,2分别存放，更新后出口时分别存放. 算法图示如下：其中过程Swap()用于交换两个参数变量的值.(iv) 计算整数步点上的坐标速度方法：procedure CalcuXv(const c:integer;const ux,uv:array of real); overload;功能：计算整数步点上的坐标速度，对应步点的系数向量由参数数组UX和UV 传入，参数c在正推时置1，逆推时置-1.算法：用(39)和(43)式，逆推时步长为-h，故可统一写为c*h.用程序语言表示如下： for j:=1 to 3 do begin p:=0;q:=0; for k:=-4 to 4 do begin p:=p+UXk+4*fj,-c*k; q:=q+UVk+4*fj,-c*k; end; p:=p/ZZFACTOR; q:=q/ZZFACTOR; xj:=p+sfj,2,3; vj:=(q+sfj,2,3 -sfj,1,3)/(c*h); end;(v) 计算中心步点上的坐标速度方法：procedure CalcuCenterXv;功能：在外推过程中校正中心步点上的坐标速度.精度较用(39)和(43)式算得的为高.算法：外推过程中sf数组中存放的是，因此要用(47)、(48)式.(vi) 计算中心步点附近非整数步点上的坐标速度方法：function CalcuXv(const nf:extended; const velosity:boolean):TVector; overload;功能：计算中心步点附近非整数步点nf上的坐标速度并返回坐标，布尔参数velosity非真时不计算速度.算法：（55）、（56）式.(iv)和(vi)的两个函数同名，但参数序列不同.这里应用了程序设计语言的函数重载功能，以便用同名函数处理不同的情况.2.3 数值积分的开头开头的任务是根据初始时刻的坐标和速度构造函数表f与和分表sf.2.3.1 开头 方法：procedure Start;功能：数值积分开头，根据初始时刻的坐标和速度构造函数表f与和分表sf.入口参数：xn个Tbody类型天体对象的数组全局变量：n天体数t0积分初始时刻h积分步长主要工作变量：sumf，sunf0迭代控制变量算法：按无摄二体运动计算所有天体在9个步点上的函数值；然后进行迭代逼近，每次迭代先逆推迭代4步，再正推迭代4步.迭代循环的控制变量sumf为序数为0的天体在-4和+4两个步点上的三个坐标之和.当sumf的相对误差小于控制量eps=10-15时结束循环，此时外推积分所需的函数值和和分表均已准备好.用Warnier图表示如下：过程调用CalcuAccel()计算每个天体的加速度.调用Iter()迭代计算函数表.2.3.2 迭代计算函数表方法：procedure Iter(const c:integer);功能：在开头过程中迭代计算函数表.积分方向取决于参数c，c=-1时逆推 (由t0向t0-h方向)，c=1 时正推(由t0向t0+h方向).程序结构：2.4 数值积分的外推外推的任务是从已构成的函数表与和分表出发，逐步外推，求出每一步点上的坐标速度，并更新函数表与和分表，直至指定时刻t.过程：procedure Sbs;功能：向前外推数值积分一步.算法：与Iter过程取参数c=1的处理相同.每步外推均先用计算的公式预估坐标速度，更新函数值表，再用计算的公式校正坐标速度，最后更新和分表.入口参数和全局变量与Start过程同.程序结构: 2.5 应用程序结构示例2.5.1 外推计算多个天体的坐标速度至指定时刻过程: procedure Evolution (tpb,pb,sw0:real);overload;功能：由t0开始，外推计算n个天体的坐标速度tpb日，外推过程中每pb日输出一次数据. 入口参数：tpb外推天数，由t0向后为正.pb输出数据间隔天数，恒为正sw0外推步长的绝对值，恒为正全局变量：n天体数t0开始外推时刻dt由t0至当前时刻的时间x天体对象的数组h外推步长,以日为时间单位.N和t0在调用过程前需要预先赋值，以下数组x各分量的初始数据也要在调用过程前预先赋值：xi.x0,xi.v0天体在t0时的坐标速度xi.mu天体的质心引力常数主要工作变量：nn由t0开始的外推步数程序结构：见下页.2.5.2 按照给定的时刻序列计算指定天体的位置方法：procedure Evolution(Table:TTable;const sw0:real;const ObNum:integer;const FieldName:string;const apparent:Boolean);overload;功能：根据数据表Table 的 T 字段给出的时刻，计算 ObNum 所给定的天体的平位置，送入序数为FieldNum 的字段,如参数 apparent 为真，同时给出视位置，送入z0字段. 参数： Table 数据表，字段为： No 序数 T 观测时刻（UT） Ra 观测赤经 Dec 观测赤纬 Z0 计算赤经 (视位置) z1 计算赤经 (真位置) O-C 观测量与计算量之差 Pem 对地球质量的偏导数 Pmm 对火星质量的偏导数 Pmx 对火星x坐标的偏导数 sw0 积分步长，恒为正 ObNum 待计算天体的序数 FieldNum 送入平位置（单位为弧度）字段的序数，apparent 视位置标志，为True时计算视位置（单位为度）送入z0字段工作变量：ObSet 需处理天体序数的集合 程序结构：见下页注意，以上两个函数也用到了函数重载功能.3 积分器的精度研究积分器精度的困难在于，对于我们所处理的多体问题，不存在分析解，无从直接比较近似数值解和严格分析解。一个可行的方法是忽略除太阳外各天体间的相互摄动，仅考虑各天体在太阳引力作用下的二体运动。此种运动存在分析解，可与数值解比较；太阳引力远大于其余各种摄动力之和，因而解的基本性质仍能保存.数值积分过程中累积误差的来源之一是2.1节提到的，每一步外推生成和分的过程中，将做为小量的右函数值的和数累加到原来的和分上去时，由于计算机字长不够，损失小量的有效位数而产生的舍入误差. 以火星为例，即使用10字节变量，坐标的字长为19位，其末位应为1018AU，当步长为0.25日时，(2)式右函数不大于105 AU,其末位应为1024AU，直接相加就可能损失5位有效数字.图1左图所示为取步长为0.25日，积分144000步（约98年），和分存储于一个10字节变量时火星的数值积分位置与分析计算位置之差.可见存在约21012 AU的累积误差. 右图所示为采用2.1节所述张家祥提出的数据结构，用两个10字节变量存储和分，取因数为220，其余条件不变时两种位置之差.可见累积误差已完全消除，只存在大小不超过1.11014的随机误差. 根据数值积分的理论，外推步数为n时，累积误差的大小约为0.1124 n3/2 2，单位同末位有效数字. 因此，外推144,000步后，随机误差约为6142003，即误差向前传递7位数字. 10字节浮点数的有效位数为十进制的19位，尚不足避免此种误差累积.引入所述数据结构，和分的有效数字增加到十进制的26位，故可以完全避免因字长不够而引起的舍入误差的积累.图1 火星的数值积分位置与分析计算位置之差. 步长为0.25日，积分144000步（约98年）左图：和分存储于一个10字节变量，有1920位有效数字右图：和分存储于两个10字节变量，因数为220，有2526位有效数字累积误差的来源之二是所用积分公式的截断误差，也即以近似公式代替精确解所产生的误差. 我们所取坐标表达式(39)中，略去了8阶和高于8阶的差分. 换言之，对于坐标，每一步积分的截断误差正比于右函数f 的8阶差分，而8阶差分又正比于8阶导数与步长8次幂之积. 选取适当的步长，即可将这一部分误差的累积减至最小. 太阳系天体中，水星所受太阳引力最大，右函数也最大，积分公式截断误差的作用相对较大. 取步长为0.25日，以两个10字节变量存储和分时，累积误差仍有21012 AU之大，如图2左图所示. 改取步长为0.125日，外推288,000步，积分相同的时间，累积误差即减小至61015 AU，如右图所示.累积误差的来源之三是计算右函数过程中因字长有限而产生的舍入误差，减小这部分误差的途径是增加积分初始数据的有效位数和程序变量的字长. 我们现在使用10字节浮点数变量和16位有效数字初值，即可将总的累积误差控制在上述范围之内.图2 水星的数值积分位置与分析计算位置之差. 和分存储于两个10字节变量，因数为220，有2526位有效数字左图：步长为0.25日，积分144000步（约98年）右图：步长为0.125日，积分288000步（约98年） 倘若考虑天体受到的所有摄动，则情况更加复杂，也没有分析解可资比较. 如上分析二体运动所得到的结论以两个10字节变量存储和分，可以基本消除更新和分时引起的舍入误差的累积；采用10字节变量计算，可以将计算过程中舍入误差引起的累积控制在可以接受的范围内对于受摄运动仍然成立. 因此需要考虑的只是与积分步长相关的积分公式截断误差引起的累积. 通过将积分步长逐步减半的方式，可以找到使总的累积误差最小的步长. 进一步减小步长时，虽然仍能减小截断误差，但由于积分步数增加，其余两种舍入误差将会占据优势，总的累积误差反会增加. 图3所示为以0.25、0.125、0.0625、0.03125和0.015625日五种步长分别积分水星坐标36000日（约98年），水星日心距和日心经度关于相邻步长的差数. 可见尽管步长成倍增加，累积误差仍然逐步减小，方法是收敛的.计算表明，水星和月球以外其他天体的坐标差，在步长由0.03125变为0.015625日时，已开始增加.下面各表详细列出了步长由0.25日开始逐步减半时，各天体位置在相邻步长之间的变化。图3 逐步缩小积分步长时，水星日心距（单位为天文单位）和日心经度（单位为弧度）在相邻步长之间的变化. ZHFACTOR=220表 步长由0.25日开始逐步减半时，各天体位置在相邻步长之间的变化ZHFACTOR=220步长 0.125日与步长0.25日的差天体日心距（AU）日心经度（mas）日心纬度（mas）水星5.510-131.210-114.310-12金星1.610-161.210-144.510-15地球4.710-115.310-112.310-11火星3.910-161.410-155.610-16木星9.910-162.310-161.410-16土星1.910-153.310-161.610-16天王星3.710-151.710-161.110-16海王星7.110-151.910-179.210-17冥王星7.310-151.410-167.810-17*月球2.810-101.810-67.910-7*太阳5.410-162.810-131.310-13地月质心1.310-155.910-157.210-15步长 0.0625日与步长0.125日的差天体日心距（AU）日心经度（mas）日心纬度（mas）水星3.410-167.410-152.810-15金星1.210-162.310-161.310-16地球2.510-152.710-151.010-15火星2.510-161.510-161.310-16木星8.910-161.310-166.510-17土星1.810-165.910-174.010-17天王星3.510-158.910-173.910-17海王星3.510-158.010-175.510-17冥王星6.710-155.710-175.810-17*月球1.210-149.210-113.710-11*太阳7.010-184.410-151.610-15地月质心2.210-161.710-169.010-17步长 0.03125日与步长0.0625日的差天体日心距（AU）日心经度（mas）日心纬度（mas）水星2.610-163.810-151.410-15金星1.210-162.110-161.210-16地球2.210-162.210-169.510-17火星2.510-161.910-161.410-16木星8.910-161.110-167.710-17土星1.810-156.910-189.810-17天王星3.510-158.510-176.110-17海王星3.510-159.610-175.510-17冥王星6.710-155.710-175.810-17*月球3.310-162.510-121.010-12*太阳2.310-172.010-148.110-15地月质心2.210-161.810-161.510-16步长 0.015625日与步长0.03125日的差天体日心距（AU）日心经度（mas）日心纬度（mas）水星1.210-161.610-155.710-16金星1.310-164.610-161.710-16地球2.310-164.210-161.510-16火星2.510-162.410-161.310-16木星9.210-161.810-161.010-16土星1.810-157.610-177.510-17天王星3.510-158.610-176.110-17海王星3.510-159.510-174.410-17冥王星6.810-157.110-175.210-17月球2.510-171.810-137.310-14太阳1.610-171.210-146.610-15地月质心2.310-164.210-161.510-16*月球为地心距、地心经度与地心纬度*太阳为太阳系质心距、太阳系质心经度与太阳系质心纬度 致谢 本文的工作得到国家天文台优秀空间天文项目基金，紫金山天文台台长基金和紫金山天文台小行星基金会等方面经费的支持. 作者衷心感谢紫金山天文台张家祥、倪维斗研究员的关心与支持，他们还详细审阅了文稿. 南大天文系学生许宏宇也阅读了文稿，并提出了好的建议，一并在此致谢.参考文献1 李广宇，倪维斗，田兰兰，PMOE 2003精密行星历表框架(I) 数学模型，紫金山天文台台刊，Vol 22,No 3-4,20032 Brouwer, D., Clemence, G. M., Methods of Celestial Mechanics, Academic Press, New York and London, 1961. 天体力学方法