第五章 小行星軌道方程計(jì)算問題——線性方程組求解的直接法

1、第五章 小行星軌道方程計(jì)算問題 線性方程組求解的直接法5.1 小行星軌道方程問題5.1.1 問題的引入一天文學(xué)家要確定一顆小行星繞太陽運(yùn)行的軌道，他在軌道平面內(nèi)建立以太陽為原點(diǎn)的直角坐標(biāo)系，其單位為天文測(cè)量單位，在5個(gè)不同的時(shí)間對(duì)小行星作了5次觀察，測(cè)得軌道上的5個(gè)點(diǎn)的坐標(biāo)數(shù)據(jù)如下表：表5.1.1 軌道上的5個(gè)點(diǎn)的坐標(biāo)數(shù)據(jù)123455.7646.2866.7597.1687.4080.6481.2021.8232.5263.360試確立小行星的軌道方程，并畫出小行星的運(yùn)動(dòng)軌線圖形。5.1.2 模型的分析 由開普勒第一定律知，小行星軌道為一橢圓，設(shè)橢圓的一般方程為： ，需要確定系數(shù)。利用已知的數(shù)

2、據(jù)，不妨設(shè)，欲確定系數(shù)等價(jià)于求解一個(gè)線性方程組： 可寫成矩陣的形式： 其中，。5.1.3 模型的假設(shè)假設(shè)：(1)小行星軌道方程滿足開普勒第一定律；(2)以上所測(cè)得數(shù)據(jù)真實(shí)有效。5.1.4 模型的建立該問題的模型為：，可見，解答上述問題就是對(duì)線性方程進(jìn)行求解。5.2 線性方程組直接解法概述直接法：利用一系列遞推公式計(jì)算有限步，能直接得到方程組的精確解。當(dāng)然，實(shí)際計(jì)算結(jié)果仍有誤差，譬如舍入誤差，舍入誤差的積累有時(shí)甚至?xí)?yán)重影響解的精度。求解線性方程組最基本的一種直接法是消去法。這是一個(gè)眾所周知的古老方法，但用在現(xiàn)代電子計(jì)算機(jī)上仍然十分有效。消去法的基本思想是，通過將一個(gè)方程乘以或除以某個(gè)常數(shù)，以及

3、將兩個(gè)方程相加減這兩種手續(xù)，逐步減少方程中的變?cè)臄?shù)目，最終使每個(gè)方程僅含一個(gè)變?cè)?，從而得出所求的解。其中高斯消去法是廣泛應(yīng)用的方法，其求解過程分為消元過程和回代過程兩個(gè)環(huán)節(jié)。消元過程將所給的方程組加工成上三角方程組。所歸結(jié)的方程組再通過回代過程得出它的解。高斯消去法由于添加了回代的過程，算法結(jié)構(gòu)稍復(fù)雜，但這種算法的改進(jìn)明顯減少了計(jì)算量。直接法比較適用于中小型方程組。對(duì)高階方程組，即使系數(shù)矩陣是稀疏的，但在運(yùn)算中很難保持稀疏性，因而有存儲(chǔ)量大，程序復(fù)雜等不足。5.3 直接解法5.3.1 高斯消去法消去法是一個(gè)古老的求解線性方程組的方法。由它改進(jìn)的選主元法是目前計(jì)算機(jī)上常用的有效的求解低階稠密矩

4、陣線性方程組的方法。例5.3.1 用消去法解方程組解 第1步，加到(5.3.2)，加到(5.3.3)，得等價(jià)方程組： 第2步，加到(5.3.5)得等價(jià)的方程組：第3步，回代法求解(5.3.6)即可求得該方程組的解為：用矩陣法描述的約化過程即為：.這種求解過程稱為具有回代的消去法。 此例可見消去法的基本思想是：用矩陣的初等行變換將系數(shù)矩陣化為具有簡單形式的矩陣（如上三角陣，單位矩陣等），而三角形方程組是很容易回代求解的。 一般的，設(shè)有個(gè)未知數(shù)的線性方程組為： （5.3.7）設(shè)則(5.3.7)化為為方便，。則消去法為：第步：計(jì)算，用乘(5.3.7)的第一方程加到第個(gè)方程中 (即實(shí)行行的初等變換)

5、。消去第2個(gè)到第個(gè)方程中的未知數(shù)得與(5.3.7)等價(jià)方程組： (5.3.8) 記為：其中(5.3.8)式中元素為進(jìn)一步需要計(jì)算的元素，公式為：。 第步，繼續(xù)上述過程消元。設(shè)第1步到第步計(jì)算已完成，得到與原方程組等價(jià)的方程組： (5.3.9)記為，下面進(jìn)行第步消元法：設(shè)，計(jì)算乘數(shù)用乘 (5.3.9)中第個(gè)方程加到第個(gè)方程消去(5.3.9)中第個(gè)方程的未知數(shù)，得到與原方程組等價(jià)的方程組： (5.3.10)記為。其中中元素計(jì)算公式為： (5.3.11)最后，重復(fù)上述過程，即且設(shè)共完成步消元計(jì)算，得到與(5.3.7)等價(jià)的三角形方程組。 (5.3.12)再用回代法求解(5.3.12)的解，計(jì)算公式為

6、： (5.3.13) 元素稱為約化的主元素。將(5.3.7)化為(5.3.12)的過程稱為消元過程。 由消元過程和回代過程求解線性方程組的方法稱為消去法。(5.3.12)的求解過程(5.3.13)稱為回代過程。定理5.3.1 （Gauss消去法）設(shè)若約化的主元素則可通過Gauss消去法（不進(jìn)行兩行的初等變換兩行交換位置）將方程組化為等價(jià)的三角形方程組(5.3.12)。消元和求解的計(jì)算公式為：消元計(jì)算 回代計(jì)算5.3.2 矩陣的三角分解下面用矩陣?yán)碚撨M(jìn)一步來分析消去法，設(shè)約化主元素由于對(duì)實(shí)行的初等變換相當(dāng)于用初等矩陣左乘，于是，消去法第1步：,有：其中：(為初等三角矩陣)消去法第步消元過程：則有

7、 (5.3.14)其中：利用遞推公式(5.3.14)則有： (5.3.15)由(5.3.15)得： (5.3.16)其中為由乘數(shù)構(gòu)成的下三角陣，為上三角矩陣，(5.3.16)表明，用矩陣?yán)碚搧矸治鱿シǎ玫揭粋€(gè)重要結(jié)果，即在條件下消去法實(shí)質(zhì)上是將分解成兩個(gè)三角矩陣的 顯然，可由消去法及行列式性質(zhì)可知，如果，則有。其中（的順序主子式）反之，可用歸納法證明：如果的順序主子式滿足：則總結(jié)以上討論，可得如下重要定理：定理5.3.2 （矩陣的三角分解）設(shè)，如果的順序主子式有，則可分解為一個(gè)單位下三角矩陣與一個(gè)上三角矩陣的乘積，即且分解是唯一的。證明 現(xiàn)僅就來證明唯一性，存在性上面已證。假若(5.3.1

8、7)且對(duì)非奇異時(shí)考慮，為單位下三角陣，為上三角陣，由假設(shè)知存在（因?yàn)榭赡婀士赡妫?，從而?5.3.17)有。上式右端為上三角陣，左端為單位下三角陣，因此左右兩端應(yīng)為單位矩陣。故即分解是唯一的。稱矩陣的三角分解（杜利特爾）分解。其中 在以上定理?xiàng)l件下，同樣可有下面的三角分解：。其中為下三角矩陣，為單位上三角矩陣，稱之為（克勞特）分解。如前例中系數(shù)矩陣的分解為：現(xiàn)設(shè)，若有分解，則(1)(2)，而求解這兩個(gè)三角形方程組是很容易的。5.3.3 Gauss消去法的計(jì)算量定理5.3.3 設(shè)為階非奇異矩陣，則用消去法解所需要的乘除法次數(shù)及加減法的次數(shù)分別為：但如果用（克萊姆）法則解，就需要計(jì)算個(gè)階行列式，若

9、行列式用子式展開，總共需要次乘法，如時(shí)消去法需要次乘除法，而克萊姆法則需要次乘法，由此可見，法則方程組的工作量太大，不便于使用。如果計(jì)算是在每秒作次乘除法計(jì)算機(jī)上進(jìn)的，那么用消去法解階方程組約需秒即可完成，而用法則大約需小時(shí)才能完成（大約相當(dāng)于年）可見，法則完全不適于在計(jì)算機(jī)上求解高維方程組。5.3.4 Gauss主元素消去法用消去法解時(shí)，設(shè)非奇異，可能出現(xiàn)，這時(shí)必須進(jìn)行行交換的Gauss消去法，但在實(shí)際計(jì)算中即使，但其絕對(duì)值很小時(shí)，用作除數(shù)，會(huì)導(dǎo)致中間結(jié)果矩陣的元素?cái)?shù)量級(jí)嚴(yán)重增長和舍入誤差的擴(kuò)散，使最后結(jié)果不可靠。例5.3.2 設(shè)有方程組解 方程組的解方法一：用消去法求解（用具有舍入的3位浮

10、點(diǎn)數(shù)進(jìn)行運(yùn)算）。 回代得解，與精確解比較，這個(gè)結(jié)果很差。方法二：用具有行交換的消去法（避免小主元）?；卮媒猓?。 這個(gè)解對(duì)于具有舍入的3位浮點(diǎn)數(shù)進(jìn)行計(jì)算，是一個(gè)很好的結(jié)果。方法一計(jì)算失敗的原因，是用了一個(gè)絕對(duì)值很小的數(shù)作除數(shù)，乘數(shù)很大，引起約化中間結(jié)果數(shù)量誤差很嚴(yán)重地增長，再舍入就使得計(jì)算結(jié)果不可靠了。這個(gè)例子告訴我們，在采用高斯消去法解方程組時(shí)，小主元可能導(dǎo)致計(jì)算失敗，故在消去法中應(yīng)避免采用絕對(duì)值很小的主元素。 對(duì)一般矩陣方程組，需要引進(jìn)主元的技巧，即在高斯消去法的每一步應(yīng)該選取系數(shù)矩陣或消元后的低階矩陣中的絕對(duì)值最大的元素作為主元素，保持乘數(shù)以便減少計(jì)算過程中的舍入誤差對(duì)計(jì)算解的影響。這個(gè)

11、例子還告訴我們，對(duì)同一個(gè)數(shù)值問題，用不同的計(jì)算方法，得到的精度大不一樣。一個(gè)計(jì)算方法，如果用此方法的計(jì)算過程中舍入誤差得到控制，對(duì)計(jì)算結(jié)果影響較小，稱此方法為數(shù)值穩(wěn)定的；否則，如果用此計(jì)算方法的計(jì)算過程中舍入誤差增長迅速，計(jì)算結(jié)果受舍入誤差影響較大，稱此方法為數(shù)值不穩(wěn)定的。因此，我們解數(shù)值問題時(shí)，應(yīng)選擇和使用數(shù)值穩(wěn)定的計(jì)算方法，否則，如果使用數(shù)值不穩(wěn)定的計(jì)算方法去解數(shù)值計(jì)算問題，就可能導(dǎo)致計(jì)算失敗。5.3.5 完全主元素消去法設(shè)有線性方程組，其中為非奇異矩陣。方程組的增廣矩陣為第1步：首先在中選主元素，即選擇使，再交換的第一行與第行元素，交換的第一列與第列元素（相當(dāng)于交換未知數(shù)與），將調(diào)到位置

12、（交換后的增廣矩陣為，其元素仍記為，然后進(jìn)行消元法計(jì)算。 第步：繼續(xù)上述過程，設(shè)已完成第步到第步計(jì)算，約化為下述形式（為簡單起見，仍記元素為元素為）：第步選主元區(qū)域 于是第步計(jì)算：對(duì)于按下述步驟從(1)計(jì)算到(3)：(1) 選主元素：選擇使；(2) 如果，則交換第行與第行元素，如果，則交換的第列與第列元素；(3) 消元計(jì)算 ；(4) 回代求解。經(jīng)過上面的過程，即從第步到第步完成選主元，交換兩行，交換兩列，消元計(jì)算。原方程組約化為：其中為未知數(shù)調(diào)換后的順序?；卮蠼猓?.3.6 列主元消去法 完全主元素消去法是解低階稠密矩陣方程組的有效方法，但完全主元素方法在選主元時(shí)要花費(fèi)一定的時(shí)間?，F(xiàn)介紹一種

13、在實(shí)際計(jì)算中常用的部分選主元（即列主元）消去法。列主元消去法，即是每次選主元時(shí)，僅依次按列選取絕對(duì)值最大的元素作為主元素，且僅交換兩行，再進(jìn)行消元計(jì)算。 設(shè)列主元消去法已經(jīng)完成第步到第步的按列選主元，交換兩行，消元計(jì)算得到與原方程組等價(jià)的方程組： 第步選主元區(qū)域第步計(jì)算如下：對(duì)于按下述步驟從(1)計(jì)算到(4)：(1) 按列主元，即確定使；(2) 如果，則為非奇異，停止計(jì)算；(3) 如果，則交換第行第行元素；(4) 消元計(jì)算： ；(5) 回代計(jì)算：計(jì)算解在常數(shù)項(xiàng)內(nèi)得到。例5.3.3 用列主元消去法求解方程組解 精確解為（舍入值）： 回代計(jì)算得到計(jì)算解：本例是具有舍入的位浮點(diǎn)數(shù)進(jìn)行運(yùn)算，所得的計(jì)算

14、解還是比較準(zhǔn)確的。 下面是完全主元素消去法框圖（圖5.3.1）：輸入 選主元素 否是否交換行輸出 否交換列且消元計(jì)算回代（當(dāng)） ）調(diào)整求解輸入計(jì)算解圖5.3.1 完全主元素消去法框圖用完全主元素消元法解，可用一整型數(shù)組開始記錄未知數(shù)次序，即，最后記錄調(diào)整后未知數(shù)的坐標(biāo)。系數(shù)陣存在二維數(shù)組內(nèi)，常數(shù)項(xiàng)存在內(nèi)，解存在數(shù)組內(nèi)。例5.3.4 若在計(jì)算過程中，只取位有效數(shù)字，試用列主元素法求解：解 第一步，選為主列元，將(5.3.18) (5.3.20)對(duì)調(diào)位置，第二步，選為列主元，不需換行， 由回代求解得：與原方程組準(zhǔn)確解比較可知，本題用位有效數(shù)字計(jì)算的列主元素法，其精度較高。 下面用矩陣運(yùn)算來描述列主

15、元素法：記 是初等排列陣（由交換單位矩陣的第行與行所得）。則列主元素法為：其中的元素滿足，由(5.3.21)得：簡記為：，其中。考察時(shí)的 (5.3.22)其中則由排列陣性質(zhì)（左乘矩陣是對(duì)矩陣進(jìn)行行變化）:已知為單位下三角陣，其元素絕對(duì)值，記。由(5.3.22)得：。其中為排列陣，為單位下三角陣，為上三角陣。這表明，對(duì)應(yīng)用列主元素法相當(dāng)于對(duì)先進(jìn)行一系列行變換后對(duì)再應(yīng)用消去法。在實(shí)際計(jì)算中我們只能在計(jì)算過程中做關(guān)于行的變換。有結(jié)論：定理5.3.4 （列主元素三角分解定理） 若為非奇異性矩陣，則存在排列矩陣使。其中為單位下三角陣，為上三角陣。 存放在的下三角部分，存放在的上三角部分。由整數(shù)型數(shù)組記錄

16、可知的情況。5.3.7 GaussJordan消去法 消去法總是消去對(duì)角線下方的元素?，F(xiàn)考慮一種修正，即消去對(duì)角線下方和上方的元素。這即為消去法。 設(shè)用消去法已完成步，于是化為等價(jià)方程組其中：在第步計(jì)算時(shí)，考慮對(duì)上述矩陣的第行上、下都進(jìn)行消元計(jì)算(1) 按列選主元素，即定義使；(2) 換行（當(dāng)）：交換第行與第行元素；(3) 計(jì)算乘數(shù) (可存放在的單元中) （）；(4) 消元計(jì)算；(5) 計(jì)算主元素。上述過程全部執(zhí)行完后有： 這表明用方法將約化為單位矩陣，計(jì)算解就在常數(shù)項(xiàng)位置得到，因此用不著回代求解。用方法解方程組的計(jì)算量大約需要次乘除法，要比消去法大些。但用方法求一個(gè)矩陣的逆矩陣還是比較合適的

17、。定理5.3.5 （法求逆矩陣）設(shè)的逆矩陣，即求階矩陣，使其中為單位矩陣，將按列寫成，為列向量，為單位列向量。于是求解，等價(jià)于求解個(gè)方程組：，所以我們可以用法求解。例5.3.5 對(duì)求。解故：5.3.8 Gauss消去法的變形(1)直接三角分解 為求解將進(jìn)行分解。即將原問題轉(zhuǎn)化為求解兩個(gè)三角形方程組：(1)，求；(2)，求。不選主元的三角分解法設(shè)，且有，其中為單位下三角陣，上三角陣。即 (5.3.23)下面說明：的元素可由步直接計(jì)算出來，其中第步定出的第行和的第列元素。由(5.3.23)有：(5.3.24)，得的第1行元素(5.3.25)，得的第1列元素到第列元素。由(5.3.23)利用矩陣乘法

18、有： 故 (5.3.26)又由(5.3.23)有：故： (5.3.27)因此可得的第行和第列的全部元素。直接分解法約需乘除法，和消去法計(jì)算量基本相當(dāng)。對(duì)計(jì)算和式，可采用雙精度累加以提高精度。選主元的三角分解法在直接三角分解中，如果計(jì)算將要中斷，或者當(dāng)絕對(duì)值很小時(shí)，按分解公式計(jì)算可能引起舍入誤差的積累。但當(dāng)為非奇異時(shí)，可通過交換的行實(shí)現(xiàn)矩陣的分解。因此可采用與列主之消去法類似的方法將直接三角分解法修改為部分選主之的三角分解法。設(shè)已完成第步分解，這時(shí)有：第步分解需要到(5.3.26) (5.3.27)兩式，為避免(5.3.27)中用小的數(shù)作除數(shù)，先引進(jìn)量：由(5.3.26)及定義，易見則由(5.3

19、.27)有：若，則我們可以用作為交換的第行與行（但我們將交換后的新元素仍記為及）于是有。控制了誤差傳播，再進(jìn)行第步分解。對(duì)一般的非奇異矩陣求逆的方法：則：。5.3.9 平方根法 利用對(duì)稱正定矩陣的三角分解得到的求解對(duì)稱正定方程組的一種有效方法平方根法。設(shè)為對(duì)稱陣，即，且的所有順序主子式均非零，則知可以唯一分解為：。為利用的非奇異性，將再分解為：為對(duì)角矩陣，為單位上三角陣，則： (5.3.30)又，由分解的唯一性即得：，代入上面(5.3.30)中得：。定理5.3.6 （對(duì)稱陣的三角分解）設(shè)為階對(duì)稱陣，且的所有順序主子式均非零。則可唯一分解為：。其中為單位下三角陣，為對(duì)角陣。 當(dāng)為對(duì)稱正定矩陣時(shí)，

20、則的所有順序主子式，而設(shè)，則于是：從而：，其中為下三角陣。定理5.3.7 （對(duì)稱正定矩陣的分解） 若為階對(duì)稱正定矩陣，則存在一個(gè)實(shí)的非奇異下三角矩陣，使。當(dāng)限定的對(duì)角元素為正時(shí)，這種分解是唯一的。 下面來考慮計(jì)算元素的公式，由，由矩陣乘法及（時(shí)）有：。于是得解對(duì)稱正定方程組的平方根法計(jì)算公式： 求解，即求解兩個(gè)三角形方程組：(1), (2)， 由(5.3.31)知，則，從而。 這表明分解過程中元素的數(shù)量級(jí)不會(huì)增長太快且對(duì)角元素恒為整數(shù)。于是不選主元素的平方根是一個(gè)數(shù)值穩(wěn)定的方法。當(dāng)求出的第列元素時(shí)，的第行亦得出，所以平方根法大約需要次乘除法，約為一般直接分解法計(jì)算量的一半。由于為對(duì)稱陣，因此在

21、計(jì)算機(jī)中只需存儲(chǔ)的下三角部分元素，共需個(gè)元素?？捎靡痪S數(shù)組存儲(chǔ)，即：。 矩陣的元素一維數(shù)組的表示為：。的元素存放在的相應(yīng)位置。 由公式(5.3.31)可見，用平方根法解對(duì)稱正定方程組時(shí)，計(jì)算的元素需要進(jìn)行開方計(jì)算，為避免開方計(jì)算，可對(duì)平方根法進(jìn)行改進(jìn)。例5.3.6 用平方根法求解：精確解為。(1)分解計(jì)算： 故，(2) 求解兩個(gè)三角方程：解，得，代入，解得：。5.3.10 追趕法在一些實(shí)際問題中常有解三對(duì)角線性方程組Ax=f的問題，即： (5.3.35)其中滿足條件：(1) (2) (3) 對(duì)于具有條件(2)的方程組(5.3.35)，我們介紹下面的追趕法求解。追趕法具有計(jì)算量少，方法簡單，算法穩(wěn)定的特點(diǎn)。定理5.3.8 設(shè)有三對(duì)角線性方程組，且滿足條件(2)，則A為非奇異矩陣。證明 用歸納法證明。顯然時(shí)，有：現(xiàn)假設(shè)定理對(duì)n-1階的滿足條件(2)的三對(duì)角矩陣成立，求證對(duì)滿足條件(2)的階三對(duì)角矩陣定理亦成立。由條件 ，則利用消元法的第步有：，顯然，其中，且有故知滿足條件(2)，利用歸納設(shè)知，故定理5.3.9 設(shè)滿足條件(2)為三對(duì)角陣。則的所有順序主子式都不為零。即：證明 由于是滿足條件(2)的階三對(duì)角陣。因此的任一個(gè)順序主子式亦滿足條件(2)的階三對(duì)角矩陣。由上一個(gè)定理即知：根據(jù)這一結(jié)論以及三角分解定理