常微分方程數值解法在子午線正反算中的應用

(武漢大學測繪學院，湖北武漢 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