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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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