C語(yǔ)言實(shí)現(xiàn)方陣的LU分解

1、對(duì) n 階方陣 A 的 LU 分解實(shí)驗(yàn)?zāi)康恼莆諏?duì) n 階方陣 A 的 LU 分解方法，并用 C 語(yǔ)言編程實(shí)現(xiàn)；觀察 L、U 矩陣的特征，并 在此基礎(chǔ)上得出求解方程方程組 AX=b 的方法。實(shí)驗(yàn)原理ij n*nij設(shè)矩陣 A = (a )的各階順序主子式均不為零，即：a(i) 0(i = 1, 2, 3L, k) ，則由定理：設(shè) n 階方陣 A 的順序主子式 D1 , D2 ,L, Dn-1 均不為零，則 A 的 LU 分解存在且唯一?？傻茫篈 有分解式 a11a12La1n a21a22La2n MMM an1 1an 2Lann u11u12Lu1 n l211 u22L u2 n = MO

2、 LU（1）OM ln1ln 2L1 unn 對(duì)比等式左邊和右邊乘積矩陣 LU 的第 r 行主對(duì)角元（含主對(duì)角元）的對(duì)應(yīng)元素，得rari = lrkukik =1(i = r,L, n; r = 1, 2,L, n).（2）再對(duì)比等式兩邊第 r 列主對(duì)角元以下（不含主對(duì)角元）的對(duì)應(yīng)元素，得rair = likukrk =1(i = r +1,L, n; r = 1, 2,L, n -1).（3）當(dāng) r=1 時(shí)，由（2）和（3）式分別解出u1i = a1i(i = 1, 2,L, n),（4）li1 =ai1u11(i = 2, 3,L, n).（5）假設(shè)已求出 U 的第 1 至 r-1 行，L

3、 的第 1 至 r-1 列，由式（2）和式（3）分別得出計(jì)算 U的第 r 行、L 的第 r 列元素的計(jì)算公式：r -1uri = ari - lrkukik =1(i = r,L, n; r = 2, 3,Ln),（6）lirr -1air - lrk ukr= k =1urr(i = r +1,L, n; r = 2, 3,L n -1).（7）稱由（4）（7）式所表示的矩陣分解為 Doolittle 分解，也稱為 LU 分解。程序?qū)崿F(xiàn)過(guò)程以下均為程序的關(guān)鍵部分：1、 輸入矩陣 Afor(i=1;iN+1;i+)%控制行數(shù)printf(請(qǐng)輸入矩陣第%d 行元素：n,i); for(j=1;j

4、N+1;j+)%控制列數(shù)scanf(%f,&Aij);%輸入每行的矩陣元素nkk2、 輸入矩陣 A 后先判斷 A 的各階順序主子式（即 k0，若不成立則不能進(jìn)行 LU 分解a ( k ) (k = 1, 2,L n ) ）是否全大于for(t=1,k=d11;tN+1;t+)%控制階數(shù)if(d11=0)&判斷矩陣第一行是否為 0，若為 0 則不分解break;eLsefor(i=t+1;iN+1;i+)%進(jìn)行第 t 次（行優(yōu)先）所有元素的消元m=dit/dtt;%計(jì)算行乘數(shù) for(j=t+1;jN;j+)dqj=dqj-m*dtj;%計(jì)算每次不為 0 的元素的值m=0;%將 m 歸零n=dt

5、+1t+1;%取主元 nk*=n;%計(jì)算各階的順序主子式3、 若可分解，進(jìn)行 LU 分解。L 為單位下三角矩陣，U 為上三角矩陣。for(r=2;rN+1;r+)%控制 U 的行數(shù)（L 的列數(shù)）for(j=r;jN+1;j+)%控制 U 的第 r 行的 j 列變化x=0;%初始化中間變量 x for(k=1;kr;k+)x+=Lrk*Ukj;%待求元素之前對(duì)應(yīng)的行和列的元素和Urj=crj-x;%計(jì)算 U 的第 r 行元素%往下是對(duì) L 的第 r 列進(jìn)行分解for(i=r+1;iN+1;i+)%控制 L 的第 r 列的 i 行變化y=0;%初始化中間變量 y for(k=1;kr;k+)y+=

6、Lik*Ukr;%待求元素之前對(duì)應(yīng)的行和列的元素和Lir=(cir-y)/Urr;%計(jì)算 L 的第 r 列元素實(shí)驗(yàn)結(jié)果及分析1、 運(yùn)行程序后，輸入待分解矩陣 A（4*4） 2100-3 -3-4-1213 A= ，得到如下結(jié)果： 123 4149-4 -132、 重新運(yùn)行程序，輸入待分解矩陣 B(4*4) 2141 3-121 B= ， 得到如下結(jié)果: 1232 5062 3、 從上述兩個(gè)不同矩陣的分解，可以看出：1） 只有矩陣為方陣且其各階順序主子式不為 0（即為非奇異矩陣），才可以進(jìn)行 LU分解，并且分解式是唯一的；否則，不能進(jìn)行 LU 分解。2） 分解的 L 矩陣為單位下三角矩陣，U 矩

7、陣為上三角矩陣。4、 實(shí)驗(yàn)簡(jiǎn)單心得：1) C 語(yǔ)言中矩陣的存儲(chǔ)同二維數(shù)組，均采用行優(yōu)先原則；本程序調(diào)用了兩個(gè)函數(shù) IsA（）和 DisA（）分別判斷 A 矩陣是否為奇異矩陣和分解 A 矩陣。2) 判斷矩陣 A 是否為奇異矩陣時(shí)，若第一個(gè)元素為 0 則其必為奇異矩陣，進(jìn)行每次 消元循環(huán)時(shí)，必須計(jì)算出每一行的行列式，而且必須設(shè)一個(gè)中間變量 n 記錄每次 消元后的主對(duì)角元素，最后通過(guò)每個(gè)主對(duì)角元素連乘 k 判斷是否為 0。3) 進(jìn)行矩陣 LU 分解時(shí)，必須同時(shí)計(jì)算 U 的第 r 行和 L 的第 r 列，在計(jì)算具體行或列 式則需先計(jì)算出 U 的第 r 行，才能計(jì)算出 L 的第 r 列。4) 本程序沒(méi)有考慮部分選主元，若考慮則必須在每次分解前判斷主元進(jìn)行行交換， 增加了程序的設(shè)計(jì)難度。5) 程序不能實(shí)現(xiàn)矩陣大小的實(shí)時(shí)改變，若要改變則必須在源程序代碼中修改宏定義 的 N 值，是一個(gè)缺陷。6) 通過(guò)程序的設(shè)計(jì)增強(qiáng)了我編程的意識(shí)，采用模塊化，若出現(xiàn)錯(cuò)誤則只需修改某個(gè) 函