[數(shù)學(xué)]數(shù)學(xué)建模中南大學(xué)數(shù)模課件

1、科學(xué)計(jì)算與數(shù)學(xué)建模中南大學(xué)數(shù)學(xué)科學(xué)與計(jì)算技術(shù)學(xué)院 小行星軌道方程計(jì)算問題2022/7/8第五章 小行星軌道方程計(jì)算問題 線性方程組求解的直接法小行星軌道方程問題 5.1線性方程組直接解法概述5.2直接解法5.3小行星軌道方程問題的模型求解5.42022/7/85.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ù)試確立小行星的軌道方程，并畫出小行星的運(yùn)動(dòng)軌線圖形。123455.7646.2866.7597

2、.1687.4080.6481.2021.8232.5263.360 5.1 小行星軌道方程問題 2022/7/85.1.2 模型的分析 由開普勒第一定律知，小行星軌道為一橢圓，設(shè)橢圓的一般方程為： ，需要確定系數(shù) 利用已知的數(shù)據(jù)，不妨設(shè) 欲確定系數(shù) 等價(jià)于求解一個(gè)線性方程組： 可寫成矩陣的形式： 2022/7/8 其中， ，2022/7/8 5.1.3 模型的假設(shè) 假設(shè)：(1)小行星軌道方程滿足開普勒第一定律；(2)以上所測(cè)得數(shù)據(jù)真實(shí)有效。 5.1.3 模型的建立 該問題的模型為： 可見，解答上述問題就是對(duì)線性方程進(jìn)行求解。 2022/7/8 5.2 線性方程組直接解法概述 直接法：利用一系

3、列遞推公式計(jì)算有限步，能直接得到方程組的精確解。當(dāng)然，實(shí)際計(jì)算結(jié)果仍有誤差，譬如舍入誤差。舍入誤差的積累有時(shí)甚至?xí)?yán)重影響解的精度 求解線性方程組最基本的一種直接法是消去法。這是一個(gè)眾所周知的古老方法，但用在現(xiàn)代電子計(jì)算機(jī)上仍然十分有效。消去法的基本思想是，通過將一個(gè)方程乘或除以某個(gè)常數(shù)，以及將兩個(gè)方程思想是，通過將一個(gè)方程乘或除以某個(gè)常數(shù)，以及將兩個(gè)方程相加減這兩種手續(xù)，逐步減少方程中的變?cè)臄?shù)目，最終使每個(gè)方程僅含一個(gè)變?cè)?，從而得出所求的解。其中高斯消去法是廣泛應(yīng)用的方法，其求解過程分為消元過程和回代過程兩個(gè)環(huán)節(jié)。消元過程將所給的方程組加工成上三角方程組。所歸結(jié)的方程組再通過回代過程得出它

4、的解。高斯消去法由于添加了回代的過程，算法結(jié)構(gòu)稍復(fù)雜，但這種算法的改進(jìn)明顯減少了計(jì)算量。 直接法比較適用于中小型方程組。對(duì)高階方程組，即使系數(shù)矩陣是稀疏的，但在運(yùn)算中很難保持稀疏性，因而有存儲(chǔ)量大，程序復(fù)雜等不足。 2022/7/85.3.1 高斯消去法 消去法是一個(gè)古老的求解線性方程組的方法。由它改進(jìn)得選主元法是目前計(jì)算機(jī)上常用的有效的求解低階稠密矩陣線性方程組的方法。例 5.3.1 用 消去法解方程組5.3 直接解法2022/7/8解：第 1 步， 加到 ， 加到 ，得等價(jià)方程組： 第 2 步， 加到 得等價(jià)的方程組：2022/7/8 第 3 步，回代法求解 即可求得該方程組的解為： 用矩

5、陣法描述的約化過程即為： 這種求解過程稱為具有回代的高斯消去法。 此例可見Gauss消去法的基本思想是：用矩陣A的初等行變換將系數(shù)矩陣化為具有簡(jiǎn)單形式的矩陣（如上三角陣，單位矩陣等），而三角形方程組是很容易回代求解的。2022/7/8 一般的，設(shè)有個(gè)未知數(shù)的線性方程組為： 2022/7/8 設(shè) 則 化為： 為方便， 則消去法為： 第1步： 計(jì)算 用 乘 的第一方程加到第 個(gè)方程中去 （即實(shí)行行的初等變換） ，消去第2個(gè)到第n個(gè)方程中的未知數(shù) 得與 等價(jià)方程組：2022/7/8 記為： 其中 式中元素 為進(jìn)一步需要計(jì)算的元素，公式為：第 ， 步，繼續(xù)上述過程消元。設(shè)第 步到第 步計(jì)算已完成，得到

6、與原方程組等價(jià)的方程組： 記為 ，下面進(jìn)行第 步消元法：2022/7/8設(shè) ，計(jì)算乘數(shù) 用 中第 個(gè)方程加到第 個(gè)方程 消去 中第 個(gè)方程 的未知數(shù) 得到與原方程組等價(jià)的方程組： 2022/7/8記為 其中 中元素計(jì)算公式為： 最后，重復(fù)上述過程，即 且設(shè) 共完成 步消元計(jì)算，得到與 等價(jià)的三角形方程組。2022/7/8 再用回代法求解 的解，計(jì)算公式為： 元素 稱為約化的主元素。將 化為 的過程稱為消元過程。 由消元過程和回代過程求解線性方程組的方法稱為Gauss 消去法。 的求解過程 稱為回代過程。2022/7/8定理（Gauss消去法）設(shè) 若約化的主元素 則可通過Gauss消去法（不進(jìn)行

7、兩行的初等變換兩行交換位置）將方程組化為等價(jià)的三角形方程組 。消元和求解的計(jì)算公式為： 1、消元計(jì)算 2、回代計(jì)算2022/7/85.3.2 矩陣的三角分解 下面用矩陣?yán)碚撨M(jìn)一步來分析 Gauss消去法，設(shè)約化主元素 由于對(duì) 實(shí)行的初等變換相當(dāng)于用初等矩陣左乘于是， 消去法第1步： , 則有： 其中： （ 為初等三角矩陣） 2022/7/8Gauss消去法第k步消元過程： 則有其中： 2022/7/8利用遞推公式則有：由 得 : 2022/7/8其中L為由乘數(shù)構(gòu)成的下三角陣,U為上三角矩陣， 表明，用矩陣?yán)碚搧矸治鯣auss消去法，得到一個(gè)重要結(jié)果，即在 條件下Gauss消去法實(shí)質(zhì)上是A將分解

8、成兩個(gè)三角矩陣的2022/7/8顯然，可由Gauss消去法及行列式性質(zhì)可知，如果 則有 其中 為順序主子式反之，可用歸納法證明：如果A的順序主子式 滿足：則 2022/7/8定理 5.3.2 （矩陣的三角分解）設(shè) ，如果 的順序主子式有 ，則 可分解為一個(gè)單位下三角矩陣與一個(gè)上三角矩陣的乘積，即 且分解是唯一的。證明 現(xiàn)僅就 來證明唯一性，存在性上面已證。假若 且對(duì) 非奇異時(shí)考慮， 為單位下三角陣， 為上三角陣，由假設(shè)知 存在（因?yàn)?可逆 故 可逆），從而由 有 ，上式右端為上三角陣，左邊為單位下三角陣，因此左右兩端應(yīng)為單位矩陣。故 即分解是唯一的。 稱矩陣的三角分解 （杜利特爾）分解。202

9、2/7/8 其中在以上定理?xiàng)l件下，同樣可有下面的三角分解： ，其中L為下三角矩陣，U為單位上三角矩陣，稱之為Crout（克勞特）分解。如前例中系數(shù)矩陣A的分解為： 現(xiàn)設(shè) ，若如分解 ，則 而求解這兩個(gè)三角形方程組是很容易的。2022/7/8定理 5.3.3 設(shè)A為n階非奇異矩陣，則用Gauss消去法解 所需要的乘除法次數(shù)及加減法的次數(shù)分別為：5.3.3 Gauss消去法的計(jì)算量但如果用 （克萊姆）法則解 ，就需要計(jì)算 個(gè)n階行列式，若行列式用子式展開，總共需要 次乘法，如 時(shí) 消去法需要430乘除法，而克萊姆法則卻需要39916800次乘法，由此可見， 法則方程組的工作量太大，不便于使用。如果

10、計(jì)算是在每秒作 次乘除法計(jì)算機(jī)上進(jìn)的，那么用 消去法解20階方程組約需0.03秒即可完成，而用 法則大約需 小時(shí)才能完成（大約相當(dāng)于 年）可見， 法則完全不適于在計(jì)算機(jī)上求解高維方程組。2022/7/8用Gauss消去法解 時(shí)，設(shè)A非奇異，但可能出現(xiàn) ，這時(shí)必須進(jìn)行行交換的Gauss消去法，但在實(shí)際計(jì)算中即使 ，但其絕對(duì)值很小時(shí)，用 作除數(shù)，會(huì)導(dǎo)致中間結(jié)果矩陣 的元素?cái)?shù)量級(jí)嚴(yán)重增長(zhǎng)和舍入誤差的擴(kuò)散，使最后結(jié)果不可靠。例5.3.2 設(shè)有方程組5.3.4 Gauss主元素消去法2022/7/8解 方程組得解方法一：用Gauss消去法求解。（用具有舍入的3位浮點(diǎn)數(shù)進(jìn)行運(yùn)算） 回代得解 ， ，與精確解

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

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

13、未知數(shù) 與 ），將 調(diào)到 位置（交換后的增廣矩陣為 ，其元素仍記為 ，然后進(jìn)行消元法計(jì)算。 5.3.5 完全主元素消去法2022/7/8第 步：繼續(xù)上述過程，設(shè)已完成第 步到第 步計(jì)算， 約化為下述形式（為簡(jiǎn)單起見，仍記 元素為 元素為 ）：域 于是第 步計(jì)算：對(duì)于 按下述步驟從 計(jì)算到第 步主元區(qū)域2022/7/8（1）選主元素：選擇 使 （2）如果 ，則交換 第 行與第 行元素，如果 ，則交換 的第 列與第 列元素。（3）消元計(jì)算2022/7/8（4）回代求解。經(jīng)過上面的過程，即從第 步到第 步完成選主元，交換兩行，交換兩列，消元計(jì)算，原方程組約化為：其中 為未知數(shù) 調(diào)換后的順序?；卮蠼?/p>

14、：2022/7/8 5.3.6 列主元消去法 完全主元素消去法是解低階稠密矩陣方程組的有效方法，但完全主元素方法在選主元時(shí)要花費(fèi)一定的時(shí)間。現(xiàn)介紹一種在實(shí)際計(jì)算中常用的部分選主元，（即列主元）消去法。列主元消去法，即是每次選主元時(shí)，僅依次按列選取絕對(duì)值最大的元素作為主元素，且僅交換兩行，再進(jìn)行消元計(jì)算。 設(shè)列主元消去法已經(jīng)完成第 步到第 步的按列選主元，交換兩行，消元計(jì)算得到與原方程組等價(jià)的方程組： 2022/7/8 第 步選主元區(qū)域2022/7/8第 步計(jì)算如下：對(duì)于 , 按下述步驟從 計(jì)算到 按列主元，即確定 使 如果 ，則 為非奇異，停止計(jì)算。 如果 ，則交換 第 行第 行元素 消元計(jì)算

15、：2022/7/8（5）回代計(jì)算： 計(jì)算解 在常數(shù)項(xiàng)內(nèi) 得到例 5.3.3 用列主元消去法求解方程組解 精確解為（舍入值）：2022/7/8回代計(jì)算得到計(jì)算解： 本例是具有舍入的位浮點(diǎn)數(shù)進(jìn)行運(yùn)算，所得的計(jì)算解還是比較準(zhǔn)確的。 下面是完全主元素消去法框圖（圖 ） 2022/7/8輸入 選主元素 否是否交換行輸出 否交換列且消元計(jì)算回代（當(dāng) ）調(diào)整求解輸入計(jì)算解圖5.3.1 完全主元素消去法框圖2022/7/8用完全主元素消元法解 ，可用一整型數(shù)組 開始記錄未知數(shù) 次序，即 ，最后記錄調(diào)整后未知數(shù)的足標(biāo)。系數(shù)陣 存在二維數(shù)組 內(nèi)，常數(shù)項(xiàng) 存在 內(nèi)，解存在數(shù)組 內(nèi)。例 5.3.4 若在計(jì)算過程中，只

16、取3位有效數(shù)字，試用列主元素法求解：2022/7/8解 第一步，選 為主列元，將 對(duì)調(diào)位置， 第二步，選 為列主元，不需換行， 2022/7/8由 回代求解得： 與原方程組準(zhǔn)確解 比較，可知，本題用3位有效數(shù)字計(jì)算的列元素法是相當(dāng)精確的。大量實(shí)踐表明：列主元法為解線性方程組的精確方法。下面用矩陣運(yùn)算來描述列主元素法記 是初等排列陣（由交換單位矩陣 的第 行與 行所得）。則列主元素法為：其中 的元素滿足 ，由 得： 2022/7/8 簡(jiǎn)記為： ，其中 。 下面考察 時(shí)的 其中， 則由排列陣性質(zhì)（左乘矩陣是對(duì)矩陣進(jìn)行行變化。）2022/7/8 已知 為單位下三角陣，其元素絕對(duì)值 ，記 。由 得：

17、。其中 為排列陣， 為單位下三角陣， 為上三角陣。這表明，對(duì) 應(yīng)用列主元素法相當(dāng)于對(duì) 先進(jìn)行一系列行變換后對(duì) 再應(yīng)用 消去法。再實(shí)際計(jì)算中我們只能在計(jì)算過程中，做關(guān)于行的變換。有結(jié)論：定理5.3.4 （列主元素三角分解定理） 若 為非奇異性矩陣，則存在排列矩陣 使 。其中 為單位下三角陣， 為上三角陣。 存放在 的下三角部分， 存放在 的上三角部分。由整數(shù)型數(shù)組記錄可知 的情況。2022/7/8 5.3.7 Gauss-Jordan 消去法 消去法總是消去對(duì)角線下方的元素?，F(xiàn)考慮一種修正，即消去對(duì)角線下方和上方的元素。這即為 消去法。 設(shè)用 消去法已完成 步，于是 化為等價(jià)方程組 其中： 20

18、22/7/8在第 步計(jì)算時(shí)，考慮對(duì)上述矩陣的第 行上、下都進(jìn)行消元計(jì)算1、按列選主元素，即定義 使2、換行（當(dāng)ikk）。交換A,b第k行與第ik行元素。3、計(jì)算乘數(shù) ( 可存放在 的單元中) （ ） 4、消元計(jì)算計(jì)算主元素 2022/7/85、計(jì)算主元素 上述過程全部執(zhí)行完后有： 這表明 用方法將 約化為單位矩陣，計(jì)算解就在常數(shù)項(xiàng)位置得到，因此用不著回代求解。用 方法接方程組的計(jì)算量大約需要 次乘除法，要比 消去法大些。但用 方法求一個(gè)矩陣的逆矩陣還是比較合適的。 2022/7/8定理 （G-J法求逆矩陣） 設(shè) 的逆矩陣 ，即求 階矩陣 ，使 其中 為單位矩陣，將 按列寫成 ， 為列向量， 為

19、單位列向量。于是求解 ，等價(jià)于求解 個(gè)方程組： ，所以我們可以用 法求解 。例 5.3.5 對(duì) 求 。2022/7/8解 2022/7/8 故： 2022/7/8 5.3.8 直接三角分解為求解 將 進(jìn)行 分解。即 將原問題轉(zhuǎn)化為求解兩個(gè)三角形方程組 求 求 。 不選主元的三角分解法設(shè) 且有 其中 為單位下三角陣， 上三角陣。即 2022/7/8下面說明： 的元素可由 步直接計(jì)算出來，其中第 步定出 的第 行和 的第 列元素。由 有： 得 的第1行元素 得 的第1列元素到第 列元素。由 利用矩陣乘法有： 故 2022/7/8 又由 有： 故： 因此可得 的第 行和第 列的全部元素。2022/7

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

21、 ）于是有 ?？刂屏苏`差傳播，再進(jìn)行第 步分解。 對(duì)一般的非奇異矩陣 求逆的方法： 則：2022/7/8 5.3.9 平方根法 利用對(duì)稱正定矩陣的三角分解而得到的求解對(duì)稱正定方程組的一種有效方法平方根法設(shè) 為對(duì)稱陣，即 ，且 的所有順序主子式均非零，則知 可以唯一分解為：為利用 的非奇異性，將 再分解為：2022/7/8 為對(duì)角矩陣， 為單位上三角陣，則： 又 ，由分解的唯一性即得： 代入上面 中得： 。 定理 5.3.6 （對(duì)稱陣的三角分解）設(shè) 為 階對(duì)稱陣，且 的所有順序主子式均非零。則 可唯一分解為： 。其中 為單位下三角陣， 為對(duì)角陣。 當(dāng) 為對(duì)稱正定矩陣時(shí)，則 的所有順序主子式 ，而

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

23、平方根法大約需要 次乘除法，約為一般直接 分解法計(jì)算量的一半。由于 為對(duì)稱陣，因此在計(jì)算機(jī)中只需存儲(chǔ) 的下三角部分元素，共需 個(gè)元素?？捎靡痪S數(shù)組存儲(chǔ)，即： 。2022/7/8 矩陣 的元素 一維數(shù)組的表示為： 。 的元素存放在 的相應(yīng)位置。 由公式 可見，用平方根法解對(duì)稱正定方程組時(shí)，計(jì)算 的元素 需要進(jìn)行開方計(jì)算，為避免開方計(jì)算，可對(duì)平方根法進(jìn)行改進(jìn)。例 5.3.6 用平方根法求解： 精確解為 。2022/7/8(1)分解計(jì)算：故，(2)求解兩個(gè)三角方程：解 得 代入 解得： 。2022/7/8 5.3.10 追趕法在一些實(shí)際問題中常有解三對(duì)角線性方程組Ax=f的問題，即：2022/7/8

24、 其中 滿足條件： 對(duì)于具有條件 的方程組 ，我們介紹下面的追趕法求解。追趕法具有計(jì)算量少，方法簡(jiǎn)單，算法穩(wěn)定的特點(diǎn)。2022/7/8定理5.3.8 設(shè)有三對(duì)角線性方程組Ax=f，且A 滿足條件 ，則A為非奇異矩陣。證明 用歸納法證明。顯然n=2時(shí)，有：現(xiàn)假設(shè)定理對(duì)n-1階的滿足條件 的三對(duì)角矩陣成立，求證對(duì)滿足條件 的 階三對(duì)角矩陣定理亦成立。由條件 ，則利用消元法的第 步有：2022/7/8 顯然 ，其中 2022/7/8且有 故知 滿足條件 ，利用歸納設(shè)知 ，故定理 5.3.9 設(shè) 滿足條件 為三對(duì)角陣。則 的所有順序主子式都不為零。即：證明 由于A是滿足條件（2）的n階三對(duì)角陣。因此A的任一個(gè)順序主子式亦滿足條件（2）的n階三對(duì)角矩陣。由上一個(gè)定理即知： 根據(jù)這一結(jié)論以及三角分解定理知，這種矩陣A可進(jìn)行三角分解：A=LU。在這里特別的有：