常微分方程数值解法在子午线正反算中的应用

(武汉大学测绘学院，湖北武汉 430000)

1 子午线正算经典算法

参考椭球具有对称性，若要求解从赤道开始到任意纬度B的子午线弧长，只需求出积分[1，3]

(1)

式中，M为子午线曲率半径，e为椭球第一偏心率，a为椭球长半轴。

为了求出M原函数，根据牛顿二项式将M展开为幂级数，然后代入式(1)中进行积分，即可得到结果。根据牛顿二项式对其进行级数展开，展至8次项得[3]

M=m0+m2sin2B+m4sin4B+m6sin6B+m8sin8B

(2)

式中

(3)

再将正弦的幂函数展开为余弦的倍数函数

(4)

将上式代入式(1)，得

M=a0-a2cos2B+a4cos4B-a6cos6B+a8cos8B

(5)

式中

将式(5)代入式(1)进行积分得

(6)

根据式(6)很容易编写出计算机程序。

2 常微分方程数值算法求解子午线弧长

(7)

子午线弧长可看做有初值的常微分方程(7)在B处的近似解。

2.1 欧拉迭代算法

对于一阶带有初值的常微分方程

(8)

在xn处，采用泰勒级数展开

y(xn+1)=y(xn+h)

略去余项，有

y(xn+1)=y(xn)+y'(xn)h

(9)

yn+1=yn+hf(xn，yn) (n=0，1，2，…)

(10)

式(10)即为欧拉公式。

2.2 欧拉-梯形迭代算法

从式(10)中不易求得yn+1，还需要在区间[xn，xn+1]上对微分方程进行积分

(11)

将式(11)右端用梯形求积公式，有

f(xn+1，y(xn+1))]

(12)

对式(12)等号右端，用近似值yn代替y(xn)，yn+1代替y(xn+1)，可得

(13)

式(13)称为梯形公式，将(10)和式(13)合用，构成如下表达式

k=0，1，2，…;n=0，1，2，…

(14)

2.3 欧拉预估-矫正算法

实际上，当h很小时，让式(14)中的梯形公式只迭代一次就结束，精度也满足要求，该式称为欧拉预估-矫正公式

k=0，1，2，…;n=0，1，2，…

(15)

2.4 龙格-库塔算法

龙格-库塔算法推导较为复杂，这里直接给出龙格-库塔算法常用的两种形式。

(1)二阶龙格-库塔算法

(16)

(2)三阶龙格-库塔算法

(17)

3 子午线弧长反算经典算法

(18)

然后开始迭代，每次都让

(19)

直到|Bi+1-Bi|<ε停止迭代，此时Bi+1即为所求的大地纬度。

4 常微分方程与数值迭代算法

根据公式(7)，可将子午线弧长与纬度看作一个带有初值的常微分方程，将数值迭代算法应用在这个常微分方程上，即可解得大地纬度B。常用的数值迭代算法有牛顿迭代、割线法以及单点迭代法，每一种迭代算法都可以与常微分方程数值解法结合使用。这里使用牛顿迭代法来进行讨论。

(20)

(21)

式(21)中的f(Bn)可由上述四种常微分数值解法求解(X已知)。因此，每次迭代都可以根据常微分方程数值解法求得每次迭代后的f(Bn)，然后进行牛顿迭代，进而求得大地纬度B。

5 程序设计与结果分析

根据上述算法，使用C#实现上述算法并设计了程序界面[8，9]，操作界面如图1所示。在此基础上实现高斯正反算及数据检验。

图1 程序主界面

通过选择不同的算法，可得到相应算法下的结果，同时，程序会给出与经典算法的差值，如图2、图3所示。

图2 子午线弧长正算算法选择

以1975国际椭球为例，分别采用上述所介绍的数值积分、常微分方程数值解法和数值迭代方法进行计算，所得结果见表1[2]。

表1 子午线弧长正算(常微分方程数值解法)

注：(1)所得子午线弧长单位均为m；(2)由于所得结果和经典算法均在米级以下，因此表格中的几种数值积分算法所得结果省去了大于km的数值。

表2 子午线弧长反算(常微分与数值迭代解算)

注：(1)所得子午线弧长单位均为m；(2)由于所得结果和经典算法只是在(")上不同，最后三列省去了度分值；(3)欧拉迭代和欧拉预估-校正公式试步长为1/1 000，二阶龙格库塔算法步长为1/100，四阶龙格库塔算法步长为1/10；(4)牛顿迭代次数为5次。

由表1可知，在子午线弧长正算中，步长1/1 000情况下的欧拉公式结果与经典算法相同，而龙格-库塔算法在迭代次数方面优于欧拉公式和经典算法。

在表2中，常微分方程数值解法所得结果与经典算法结果基本一致，最大相差0.006 7″(基本可以忽略此差值)，并且牛顿迭代法与欧拉迭代算法相结合，弥补了欧拉公式精度低且步长小的缺点。

6 结束语

首先验证了数值积分[1]和数值迭代[2]算法在子午线正反算中的正确性，并在此基础上使用欧拉迭代、欧拉预估-矫正、龙格库塔三种常见的常微分数值解法对子午线弧长进行正反算，并与传统的子午线正反算结果进行比较。

基于公式推导及计算结果，常微分数值解法结果和传统算法结果基本一致，并且具有实现简单，迭代次数少、速度快等优点。

[1] 郑红晓，张红方，雷伟伟.子午线弧长计算的数值积分算法及其比较[J].铁道勘察，2014，40(6):8-10

[2] 郑红晓，张红方，雷伟伟.计算底点纬度Bf的数值迭代算法及其比较[J].测绘与空间地理信息，2015，38(2):42-44

[3] 孔祥元，郭际明，刘宗泉.大地测量学基础[M].武汉：武汉大学出版社，2006

[4] 严伯铎.椭球子午线弧长的一种计算方法[J].地矿测绘，2003(3):7-10

[5] 李信真，车刚明，欧阳洁，等.计算方法[M].西安：西北工业大学出版社，2010

[6] 利庆扬，王能超，毅大义，等.数值分析[M].北京:清华大学出版社，2001

[7] 严伯铎.椭球子午线弧长的一种计算方法[J].地矿测绘，2003(3):7-10

[8] JonSkeet.深入理解C#[M].姚琪琳，译.北京：人民邮电出版社，2014

[9] 里克特.CLR via C#[M].周靖，译.北京：清华大学出版社，2010

[10] 易维勇，边少锋，朱汉泉.子午线弧长的解析型幂级数确定[J].测绘学院学报，2000(3):167-171

[11] 牛卓立.以空间直角坐标为参数的子午线弧长计算公式[J].测绘通报， 2001(11):14-15

[12] 过家春.子午线弧长公式的简化及其泰勒级数解释[J].测绘学报，2014，43(2):125-130