子午线弧长的递归算法

高新疆

摘要：本文在分析子午线弧长公式结构的基础上，提出了子午线弧长的递归算法，并设计了实现算法的计算机程序。通过与数值积分算法的比较，证明了递归算法无论是计算速度、计算精度，还是计算结果的稳定性，都要优于数值积分算法。

主题词：子午线弧长，递归算法

一、引言

在大地测量中，子午线弧长的计算是一项基本的工序，其计算结果的精度决定了坐标正反算、梯形图幅边长计算的精度。

当计算从赤道至点P的子午线弧长时，公式如下[1]：

式中，为参考椭球的长半轴

为参考椭球的第一偏心率

B为P点的纬度

由子午弧长公式看出，子午线弧长的计算精度取决于积分式的计算精度。一般情况下，该积分式都是选用多项式逼近积分式进行积分计算，称为数值积分算法。常用的数值积分算法有梯形算法、抛物线算法（Simpson公式）、Romberg方法（逐次分半加速法）及其复化算法等。

但通过对子午线弧长公式内的积分式进行分析，可以发现对于这个特定的积分式，可以用另一种方法进行计算，且效果比数值积分要好。具体分析及结论如下：

二、对直接积分法的数值分析

积分因子按二项展开式展开，得[1]：

上式也符合递归算法的形式，也可用递归算法来实现。

需要说明的是，由于子午线弧长是由积分式与一个相对较大的数相乘而得，这个较大的数相当于一个放大器，可以将积分式计算中较小的误差放大。比如说，若积分式计算的精度为10-10，则计算出的子午线弧长精度只能达到微米。为了计算出高精度的子午线弧长，必须提高积分式的计算精度。

三、子午线弧长积分递归算法的程序设计

如上所述，子午线弧长公式积分式由三部分相乘而得，可分别对（3）式中的bi、Si两部分进行编程计算，而第一偏心率的幂，由于计算简单，可不必单独编程。

1、二项式系数的计算：

二项式系数的编程计算可用递归来进行，程序源代码如下：

longdoubleBinomial(longdoublepower,intterm)

{

longdoublebinResult=1;

if(term==1)

binResult=power;

else

binResult=(power-term+1)*Binomial(power,term-1)/term;

returnbinResult;

}

程序中，参数power为二项展开式中的幂，即式（3）中的-3/2；参数term为二项式系数的项序，即式（3）中的i。

2、正弦函数幂的积分运算：

正弦函数幂的积分递推公式为[3]：

以上公式也可用递归算法来实现，程序源代码如下：

longdoubleIntegrateSinN(longdoubleradAngle,intpower)

{

longdoubleintegResult;

if(power==0)integResult=radAngle;

elseif(power==1)integResult=1-cos(radAngle);

elseintegResult=(power-1)*IntegrateSinN(radAngle,power-2)/power-pow(sin(radAngle),power-1)*cos(radAngle)/power;

returnintegResult;

}

程序中，参数radAngle为以弧度为单位的角度值，在这里指的是纬度；参数power为正弦函数的幂次，也即级数项序的2倍。

3、子午弧长的计算(取项至8，即16次项)

longdoubleMeridian(longdoublelatitude,longdoublelongradiu,longdoublesquare_e)

{

longdoublemerResult;

intcounter=1;

merResult=IntegrateSinN(latitude,0);

for(counter=1;counter<8;counter++)merResult+=pow(-1,counter)*Binomial(-1.5,counter)*pow(square_e,counter)*IntegrateSinN(latitude,counter*2);

merResult*=longradiu*(1-square_e);

returnmerResult;

}

程序中，参数longradiu为椭球的长半轴，参数square_e为椭球第一偏心率的平方。该程序也可用两次迭代结果的差值来判断程序是否结束，限于篇幅，不再列出源代码。

以上程序，在WindowsXP系统下、VisiualC++6.0环境下运行通过。

四、数值积分法与直接积分法的算法比较

为了比较直接积分法与数值积分法的优劣，笔者编写了复化simpson公式[2]数值积分和Romberg方法（逐次分半加速法）[2]数值积分程序，分别用递归算法、复化simpson公式积分法及Romberg数值积分法计算了纬度40°的子午线弧长及子午线弧长公式中的积分式，并统计了计算用时，计算结果详见表1、表2及表3，参考椭球为1975年国际椭球。

表1递归算法子午线弧长计算表

由上表可以看出，递归算法收敛较快，只需计算到第7项，即14次幂的项，就可使积分式达到10-16的精度，而此时子午线弧长的精度也达到了10-10。再以增加项数来提高计算精度，一是受到计算机浮点数精度的限制而很难实现，同时也没有现实意义。

表2复化simpson公式积分法子午线弧长计算表

表中区间数用以控制计算的精度，对于复化simpson公式，实际分为2n个区间。即：步长=（积分上限-积分下限）/2n。

从表中所列数据可以看出，用复化simpson公式积分法计算收敛很慢，计算精度达到10-2时，计算用时为13分钟。因此说这种方法的计算精度不高。

表3Romberg数值积分法子午线弧长计算表

由表中数据可以看出，Romberg数值积分法的计算速度要比复化simpson公式快得多，与递归算法的计算速度相当。但计算的稳定性要劣于递归算法。积分式小数点后第15位之后的数值存在着摆动现象。也就是说，Romberg数值积分法计算子午线弧长只能精确到小数点后第7位，也即计算精度只能到10-7。

五、结论

从以上的分析中可以得出如下结论：

1、在相同的计算时间内，递归算法计算子午线弧长的精度明显要比数值积分的精度高；

2、递归算法的收敛速度要比数值积分快，也即若计算结果精度相同，递归算法的计算用时要少于数值积分算法。

3、递归算法的稳定性要优于数值积分法。

参考文献

[1]同济大学主编．高等数学（上、下册）[M]．北京：高等教育出版社1980，第二版．

[2]邓建中、葛仁杰、程正兴．计算方法[M]．西安：西安交通大学出版社，1985．

作者简介：

高新疆，汉族，1975年12月生，工程师，籍贯:山东省栖霞县，现从事工程测量、数字化测图工作。身份证号：650108197512231435。

联系地址：新疆维吾尔自治区昌吉市健康东路9号

邮编：831100电话：[phone]