常微分方程數(shù)值解法在子午線正反算中的應(yīng)用

(武漢大學(xué)測(cè)繪學(xué)院，湖北武漢 430000)

1 子午線正算經(jīng)典算法

參考橢球具有對(duì)稱性，若要求解從赤道開始到任意緯度B的子午線弧長(zhǎng)，只需求出積分[1，3]

(1)

式中，M為子午線曲率半徑，e為橢球第一偏心率，a為橢球長(zhǎng)半軸。

為了求出M原函數(shù)，根據(jù)牛頓二項(xiàng)式將M展開為冪級(jí)數(shù)，然后代入式(1)中進(jìn)行積分，即可得到結(jié)果。根據(jù)牛頓二項(xiàng)式對(duì)其進(jìn)行級(jí)數(shù)展開，展至8次項(xiàng)得[3]

M=m0+m2sin2B+m4sin4B+m6sin6B+m8sin8B

(2)

式中

(3)

再將正弦的冪函數(shù)展開為余弦的倍數(shù)函數(shù)

(4)

將上式代入式(1)，得

M=a0-a2cos2B+a4cos4B-a6cos6B+a8cos8B

(5)

式中

將式(5)代入式(1)進(jìn)行積分得

(6)

根據(jù)式(6)很容易編寫出計(jì)算機(jī)程序。

2 常微分方程數(shù)值算法求解子午線弧長(zhǎng)

(7)

子午線弧長(zhǎng)可看做有初值的常微分方程(7)在B處的近似解。

2.1 歐拉迭代算法

對(duì)于一階帶有初值的常微分方程

(8)

在xn處，采用泰勒級(jí)數(shù)展開

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

略去余項(xiàng)，有

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，還需要在區(qū)間[xn，xn+1]上對(duì)微分方程進(jìn)行積分

(11)

將式(11)右端用梯形求積公式，有

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

(12)

對(duì)式(12)等號(hào)右端，用近似值yn代替y(xn)，yn+1代替y(xn+1)，可得

(13)

式(13)稱為梯形公式，將(10)和式(13)合用，構(gòu)成如下表達(dá)式

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

(14)

2.3 歐拉預(yù)估-矯正算法

實(shí)際上，當(dāng)h很小時(shí)，讓式(14)中的梯形公式只迭代一次就結(jié)束，精度也滿足要求，該式稱為歐拉預(yù)估-矯正公式

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

(15)

2.4 龍格-庫(kù)塔算法

龍格-庫(kù)塔算法推導(dǎo)較為復(fù)雜，這里直接給出龍格-庫(kù)塔算法常用的兩種形式。

(1)二階龍格-庫(kù)塔算法

(16)

(2)三階龍格-庫(kù)塔算法

(17)

3 子午線弧長(zhǎng)反算經(jīng)典算法

(18)

然后開始迭代，每次都讓

(19)

直到|Bi+1-Bi|<ε停止迭代，此時(shí)Bi+1即為所求的大地緯度。

4 常微分方程與數(shù)值迭代算法

根據(jù)公式(7)，可將子午線弧長(zhǎng)與緯度看作一個(gè)帶有初值的常微分方程，將數(shù)值迭代算法應(yīng)用在這個(gè)常微分方程上，即可解得大地緯度B。常用的數(shù)值迭代算法有牛頓迭代、割線法以及單點(diǎn)迭代法，每一種迭代算法都可以與常微分方程數(shù)值解法結(jié)合使用。這里使用牛頓迭代法來(lái)進(jìn)行討論。

(20)

(21)

式(21)中的f(Bn)可由上述四種常微分?jǐn)?shù)值解法求解(X已知)。因此，每次迭代都可以根據(jù)常微分方程數(shù)值解法求得每次迭代后的f(Bn)，然后進(jìn)行牛頓迭代，進(jìn)而求得大地緯度B。

5 程序設(shè)計(jì)與結(jié)果分析

根據(jù)上述算法，使用C#實(shí)現(xiàn)上述算法并設(shè)計(jì)了程序界面[8，9]，操作界面如圖1所示。在此基礎(chǔ)上實(shí)現(xiàn)高斯正反算及數(shù)據(jù)檢驗(yàn)。

圖1 程序主界面

通過(guò)選擇不同的算法，可得到相應(yīng)算法下的結(jié)果，同時(shí)，程序會(huì)給出與經(jīng)典算法的差值，如圖2、圖3所示。

圖2 子午線弧長(zhǎng)正算算法選擇

以1975國(guó)際橢球?yàn)槔?，分別采用上述所介紹的數(shù)值積分、常微分方程數(shù)值解法和數(shù)值迭代方法進(jìn)行計(jì)算，所得結(jié)果見表1[2]。

表1 子午線弧長(zhǎng)正算(常微分方程數(shù)值解法)

注：(1)所得子午線弧長(zhǎng)單位均為m；(2)由于所得結(jié)果和經(jīng)典算法均在米級(jí)以下，因此表格中的幾種數(shù)值積分算法所得結(jié)果省去了大于km的數(shù)值。

表2 子午線弧長(zhǎng)反算(常微分與數(shù)值迭代解算)

注：(1)所得子午線弧長(zhǎng)單位均為m；(2)由于所得結(jié)果和經(jīng)典算法只是在(")上不同，最后三列省去了度分值；(3)歐拉迭代和歐拉預(yù)估-校正公式試步長(zhǎng)為1/1 000，二階龍格庫(kù)塔算法步長(zhǎng)為1/100，四階龍格庫(kù)塔算法步長(zhǎng)為1/10；(4)牛頓迭代次數(shù)為5次。

由表1可知，在子午線弧長(zhǎng)正算中，步長(zhǎng)1/1 000情況下的歐拉公式結(jié)果與經(jīng)典算法相同，而龍格-庫(kù)塔算法在迭代次數(shù)方面優(yōu)于歐拉公式和經(jīng)典算法。

在表2中，常微分方程數(shù)值解法所得結(jié)果與經(jīng)典算法結(jié)果基本一致，最大相差0.006 7″(基本可以忽略此差值)，并且牛頓迭代法與歐拉迭代算法相結(jié)合，彌補(bǔ)了歐拉公式精度低且步長(zhǎng)小的缺點(diǎn)。

6 結(jié)束語(yǔ)

首先驗(yàn)證了數(shù)值積分[1]和數(shù)值迭代[2]算法在子午線正反算中的正確性，并在此基礎(chǔ)上使用歐拉迭代、歐拉預(yù)估-矯正、龍格庫(kù)塔三種常見的常微分?jǐn)?shù)值解法對(duì)子午線弧長(zhǎng)進(jìn)行正反算，并與傳統(tǒng)的子午線正反算結(jié)果進(jìn)行比較。

基于公式推導(dǎo)及計(jì)算結(jié)果，常微分?jǐn)?shù)值解法結(jié)果和傳統(tǒng)算法結(jié)果基本一致，并且具有實(shí)現(xiàn)簡(jiǎn)單，迭代次數(shù)少、速度快等優(yōu)點(diǎn)。

[1] 鄭紅曉，張紅方，雷偉偉.子午線弧長(zhǎng)計(jì)算的數(shù)值積分算法及其比較[J].鐵道勘察，2014，40(6):8-10

[2] 鄭紅曉，張紅方，雷偉偉.計(jì)算底點(diǎn)緯度Bf的數(shù)值迭代算法及其比較[J].測(cè)繪與空間地理信息，2015，38(2):42-44

[3] 孔祥元，郭際明，劉宗泉.大地測(cè)量學(xué)基礎(chǔ)[M].武漢：武漢大學(xué)出版社，2006

[4] 嚴(yán)伯鐸.橢球子午線弧長(zhǎng)的一種計(jì)算方法[J].地礦測(cè)繪，2003(3):7-10

[5] 李信真，車剛明，歐陽(yáng)潔，等.計(jì)算方法[M].西安：西北工業(yè)大學(xué)出版社，2010

[6] 利慶揚(yáng)，王能超，毅大義，等.數(shù)值分析[M].北京:清華大學(xué)出版社，2001

[7] 嚴(yán)伯鐸.橢球子午線弧長(zhǎng)的一種計(jì)算方法[J].地礦測(cè)繪，2003(3):7-10

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

[9] 里克特.CLR via C#[M].周靖，譯.北京：清華大學(xué)出版社，2010

[10] 易維勇，邊少鋒，朱漢泉.子午線弧長(zhǎng)的解析型冪級(jí)數(shù)確定[J].測(cè)繪學(xué)院學(xué)報(bào)，2000(3):167-171

[11] 牛卓立.以空間直角坐標(biāo)為參數(shù)的子午線弧長(zhǎng)計(jì)算公式[J].測(cè)繪通報(bào)， 2001(11):14-15

[12] 過(guò)家春.子午線弧長(zhǎng)公式的簡(jiǎn)化及其泰勒級(jí)數(shù)解釋[J].測(cè)繪學(xué)報(bào)，2014，43(2):125-130