【正文】 0 0 計算方法實驗報告 7 結(jié)果分析 輸入 A 非奇矩陣， 與計算結(jié)果相同，達到了要求， 可以實現(xiàn) 對 N 階非奇矩陣的 LU 分解。 計算方法實驗報告 8 問題 3 三、 編寫程序?qū)崿F(xiàn)大規(guī)模方程組的列主元高斯消去法程序，并對所附的方程組進行求解。具體的要求參見所附的相關(guān)文檔。 針對本專業(yè)中所碰到的實際問題，提煉一個使用方程組進行求解的例子，并對求解過程進行分析、求解。 問題分析 有線性方程組 Ax b? ，設(shè) A 是可逆矩陣。 列主元消去 法的基本思想就是通過列主元的選取將初等行變換作用于方程組的增廣矩陣 ? ?|B A b? ，將其中的 A 變換成一個上三角矩陣，然后求解這個三角形方程組。 將方程組用增廣矩陣 ? ? ? ? ( 1 )| ij nnB A b a ???? 表示 。 第一步 ：消元過程，對 1, 2, , 1kn?? 選主元，找 ? ?, 1, ,ki k k n?? 使得 , maxki k ikk i naa??? 如果 , 0kika ? ，則矩陣 A 奇異，程序結(jié)束；否則執(zhí)行（ 3）； 如果 kik? ，則交換第 k 行與第 ki 行對應(yīng)元素位置，kkj i jaa?， , , 1j k n??； 消元，對 ,i k n? ，計算 /,ik ik kkl a a? 對 1, , 1j k n? ? ?，計算 .ij ij ik kja a l a?? 第二步 ：回代過程： 若 0,nna ? 則矩陣奇異，程序結(jié)束；否則執(zhí)行（ 2）； ,1/。n n n nnx a a?? 對 1, , 2,1in?? ，計算 計算方法實驗報告 9 ,1 1 /ni i n ij j iijix a a x a? ??????????? 對較小規(guī)模的系數(shù)矩陣 A 可按正常方法求解，但對于 所給的 43200 階矩陣，不可能有那么大的儲存空間來儲存，考慮到帶狀矩陣的算法不對零元進行計算，因此只需要儲存非零元的即可?？刹捎?“ 對角線 ” 存儲方式，按下列方式存放。其中下帶寬 p=1，上帶寬 q=3. ( 0 a11 a12 a13 a14a21 a22 a23 a24 a25a32 a33 a34 a35 a36? ? ? ? ?an,n?1 ann 0 0 0 ) 它需要 n （ p +q +1） 個數(shù)組，節(jié)省了很大存儲空間，適合解決大規(guī)模帶狀壓縮矩陣的問題， fun00 fun005 就是按照這種方式求解的。按這種方式存放，矩陣的 （ i， j）元在數(shù)組中的（ i， ji+p+1）位置。 Matlab 求解程序見附錄三。 實驗結(jié)果 調(diào)試程序用 id =F1E1D1A0 ver =102 n = 10 方程組的解： x[1]= x[2]= x[3]= x[4]= x[5]= x[6]= x[7]= x[8]= 計算方法實驗報告 10 x[9]= x[10]= 用時 50s。 id =F1E1D1A0 ver =101 n = 1024 計算方法實驗報告 11 用時 75s。 id =F1E1D1A0 ver =102 n =1560 測試矩陣 id =F1E1D1A0 ver =201 n =10 計算方法實驗報告 12 p =1 q=3 方程組的解： x[1]= x[2]= x[3]= x[4]= x[5]= x[6]= x[7]= x[8]= x[9]= x[10]= 將系數(shù)矩陣轉(zhuǎn)換成正常矩陣存儲，驗證結(jié)果一樣。 用時 60s。 id =F1E1D1A0 ver =201 n =43200 計算方法實驗報告 13 p =1 q =3 結(jié)果分析 由上面結(jié)果可以看出，所編寫程序可以解決正常系數(shù)矩陣方程組，同時能解決帶狀系數(shù)矩陣方程組，并且在解決該類系數(shù)矩陣時具有很高的效率。 工程實例 在工程技術(shù)中所遇到的電路，大多數(shù)是很復(fù)雜的，這些電路是由電器元件按照一定方式互相連接而構(gòu)成的網(wǎng)絡(luò)。在電路中，含有元件的導線稱為支路，而三條或三條以上的支路的會合點稱為節(jié)點。以下圖為例，求解電路網(wǎng)絡(luò)中各條支路計算方法實驗報告 14 上的電流。 解 : a) 問題分析： 對于這類問題的計算，通常采用基爾霍夫（ Kirchhoff）定律來解決。設(shè)各節(jié)點的電流如圖所示，則由基爾霍夫第一定律（簡記為 KCL） 回路電壓的 關(guān)系 ，任何一個閉合回路的電流服從希爾霍夫電壓定律：沿某個方向環(huán)繞回路一周的所有電壓降 U 的代數(shù)和等于沿同一方向環(huán)繞該回路一周的電源電壓的代數(shù)和。 b) 算法： 回路 1 的電路方程為 12I1 ?4I2 ?7I3 = 40。 回路 2 的電路方程為 ?4I1 +13I2 ?5I4 = ?10。 回路 3 的電路方程為 ?7I1 +15I3 ?6I4 = 30。 回路 4 的電路方程為 ?5I2 ?6I3 + 14I4 = 20. 于是求各個支路的電流就歸結(jié)為下面齊次線性方程組的求解 {12I1 ?4I2 ? 7I3 = 40?4I1 +13I2 ? 5I4 = ?10?7I1 + 15I3 ?6I4 = 30?5I2 ?6I3 +14I4 = 20 c) Matlab 求解 ： 利用列主元高斯消去法函數(shù) clc。 clear all。 A=[12,4,7,0。4,13,0,5。7,0,15,6。0,5,6,14]。 b=[40。10。30。20]。 4 0 V3 0 V2 0 V1 0 V4 Ω5 Ω7 Ω1 Ω 4 Ω6 Ω 3 Ω2 ΩI1I2I3I4ABCDE計算方法實驗報告 15 gauss(A,b) 方程組的解： x[1]= x[2]= x[3]= x[4]= 計算方法實驗報告 16 問題 4 四、 已知某產(chǎn)品從 1900 年到 2020 年每隔 10 年的產(chǎn)量 ,用多項式插值和三次樣條插值的方法，畫出每隔一年的插值曲線的圖形 , 試計算并比較在不同方法下的 2020 年以及 2020 年的產(chǎn)量。 年份 1900 1910 1920 1930 1940 1950 產(chǎn)量 年份 1960 1970 1980 1990 2020 2020 產(chǎn)量 問題分析 a) Newton 多項式 插值 法 Newton 插值多項式的表達式如下： 0 1 0 0 1 1( ) ( ) ( ) ( ) ( )n n nN x c c x x c x x x x x x ?? ? ? ? ??? ? ? ? ??? ? 其中每一項的系數(shù) ci的表達式如下： 1 2 0 1 1010[ , , , ] [ , , , ][ , , , ] iiiiif x x x f x x xc f x x x xx ???? ? ???? ??? ? ? 即為 f (x)在點 01, , , ix x x??? 處的 i 階差商，（ ? ? ()iif x f x? ， 1, 2 , ,in? ），由差商 01[ , , , ]if x x x??? 的性質(zhì)可知： ? ?01 0 0 1[ , , , ] ( ) iiijj k jkkjf x x x f x xx? ????? ? ?? ? 計算方法實驗報告 17 牛頓插值的程序?qū)崿F(xiàn)方法： 第一步：計算 ? ? ? ? ? ? ? ?0 0 1 0 1 2 0 1 2, , , , , , , nf x f x x f x x x f x x x x、 、 、 、。 第二步：計算牛頓插值多項式中 01[ , , , ]if x x x??? 0 1 1( ) ( ) ( )ix x x x x x?? ? ??? ?，1, 2 , ,in? ，得到 n 個多項式。 第三步：將第二步得到的 n 個多項式相加，得到牛頓插值多項式。 第四步：利用所得到的插值多項式，估算 x 取其它值時 ??fx的值。 第五步：作出所求多項式在插值結(jié)點周圍的函數(shù)圖像。 b) 三次樣條函數(shù)插值法 則其必滿足的三次樣條插值函數(shù)是如果 ,)()( xfxS 插值條件： njyxS jj ,1,0,)( ??? 連續(xù)性條件： 1,1,)()(l i m ????? njyxSxS jjxx j ? 一階導數(shù)連續(xù)條件： 1,1,)()(l i m ??????? njmxSxS jjxx j ? 二階導數(shù)連續(xù)條件： 1,1),()(l i m ???????? njxSxS jxx j ? 因為 s(x)在每個小區(qū)間上是一個次小于三次的多項式，故有四個未知系數(shù)；因為 s(x)有 n 分段，從而共有 4n 個未知系數(shù)，不僅計算量大，而且還占用很多計算機存儲，因此用下面方法求解 )(xS 。 選擇二階導數(shù)作為待定參數(shù) : 由于三次樣條 S(x)是三次多項式 ,故它的二階導數(shù)是一次多項式 ,從而 ( ) ( 0 , 1 , 2 , , )kkM S x k n????( 1 ) 1111() ( ) ( ) , [ , ]()k k k kkk k k k kk k k kS x M MMS x x x M x x xS x M h?????? ?? ???? ? ? ? ?? ?? ??1 ( 0 , 1 , , 1 )k k kh x x k n?? ? ? ?( 2 ) 21( ) ( ) ( )2kkk k k k kkMMS x x x M x x ph? ??? ? ? ? ? ?( 3 ) 321( ) ( ) ( ) ( )62k k kk k k k k kkM M MS x x x x x p x x qh? ?? ? ? ? ? ? ? ?計算方法實驗報告 18 這個基本方程只包括了 n1 個方程 !但卻有 n 個二階導數(shù)需要待定 ,這是一個欠定方程組 ,還需要根據(jù)邊界條件再確定兩個方程 i. 自然樣條 ii. 已知兩端點的一階導數(shù)值 ()k k k k kS x y q y? ? ?11 1 1 1 1 2( ) ( ) [ , ] 6kkk k k k k k k k kMMS x S x y p f x x h?? ? ? ? ? ?? ? ? ? ?321 2 3 4( ) ( ) ( ) ( ) ( 2 .7 )k k k k k k kS x s x x s x x s x x s? ? ? ? ? ? ?1121314622[ , ]6kkkkkkkkk k k k kk k kMMshMsMMs p f x x hs q y????????????? ?? ? ? ??? ???0 0, 0nMM?? 1 1 12 2 2 22 2 2 21 1 126262626n n n nn n nMdMdMdMd??????? ? ? ?? ? ?? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ??? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ?00( ) , ( )nnm f x m f x????0 0 0 0 1 0 0 01 1 1( ) 2 6 ( ) /( ) 2 6 ( ) /n n n n n n n nS x m M M d m hS x m M M m d h? ? ?? ? ? ? ??????? ? ? ? ???計算方法實驗報告 19 當然還有其他的形式 。 三次樣條插值一般情況只能評估插值點區(qū)域內(nèi)對求函數(shù)值，本算法采用第一類邊界條件的情況下，對最后兩個插值點進行延伸， 即對 20202020 間插值的時候采用 20202020 間的 M 值，計算出插值函數(shù)， 來求 2020 年的函數(shù)值。 Matlab 求解程序見附錄四。 實驗結(jié)果 a) Newton 多項式 插值 法實驗結(jié)果 Newton 插值 2020 年的產(chǎn)量及插值函數(shù) Nx = (15977*x)/10000(1119*(x1900)*(x 1910))/100000 + (4604768489399731*(x 1900)*(x1910)*(x1920))/46116