矩陣LU分解求逆詳細分析與C語言實現(xiàn)

1、題目要求給定一個多維矩陣，實現(xiàn)該矩陣的求逆運算。1、理論分析矩陣的一種有效而廣泛應(yīng)用的分解方法是矩陣的 LU三角分解， 將一個n階矩陣A分解為一個下三角矩陣L和一個上三角矩陣U的乘 積。所以首先對矩陣進行三角分解，這里采用 Doolittle分解，即分解 為一個下三角矩陣（對角元素為1）,和一個上三角矩陣的乘積。再進 行相應(yīng)的處理。所以，矩陣求逆的算法流程可表述如下：'輸輸A輸輸輸輸|圖1矩陣求逆流程圖1）進行LU分解;2）對分解后的L陣（下三角矩陣）和U陣（上三角矩陣）進行求逆；;3）L陣的逆矩陣和U陣的逆矩陣相乘，即可求得原來矩陣的逆。即：A -1 = （ LU ） -1 = U

2、-1 L-1（1）1.1矩陣的LU分解若n階方陣A"CnCn*n的各階順序主子式不等于零，即:a11 a12-a1ka21 a22-a2kk =：:i'a#0,(k= 1,2，,n),ak1 ak2-akk則A的LU分解A=LxU存在且。9ar.-an9A = ar1 9ar.-arn9an1 a nr* -a nn _-19+ .1U 11 Un Uml+ .9Lr1 1aa+U U=rr= rn+ .a=LUJ LnrW1U nn -由矩陣的乘法原理,可推導(dǎo)出LU分解的迭代算法U0j = a0j,( j=0,1,2，,n-1),(2)(3)(4)18(5)(6)a”lj&

3、#176; =工,(| =0,1,2, ,n-1),U00r -1arj lik Ukjr -1k =1 - ad- IrkUj , lrrk T(r =0,1,2, ,n-1;j =r, ,n-1),r -1ar- likUkr|. = 1 ir,Urr(r =0,1,2, ,n -1;i =r 1, ,n -1)矩陣的LU分解是一個循環(huán)迭代的過程，U矩陣是從第1行迭代到第n 行,而L矩陣則是從第1列迭代到第n列，且U矩陣先于L矩陣一個 節(jié)拍。1.2 L矩陣和U矩陣求逆首先假設(shè)下三角矩陣L的逆矩陣為l ,不失一般性，考慮4階的情況,利用Ll =1 ,有:(1) look 1 = L 1 ,1

4、 1 22 = L22 ,1 33 = L33 ；(2) "0 = T00 (|11L10 )；(4)(3) |2 = _|00(|21L1 l22 L20)；30 = r00 (31、0 +32 L20 *r33 L30)。從而求得下三角矩陣L的逆矩陣R式如下:j廣j(8)七=-"叭)，i j k=i 1.0,ij上三角矩陣U的逆矩陣可以由下式得到:ujlI(9)('UjkUki),Ijk=j 1.，l j0矩陣求逆是一個迭代的過程，依次循環(huán)，迭代n1次，求出整個逆矩陣。其中U矩陣的循環(huán)迭代時按行順序，列倒序進行， L矩陣的 循環(huán)迭代按列順序，行順序進行，直到計算

5、出整個矩陣的所有結(jié)果為 止。1.3矩陣相乘上三角矩陣U的逆矩陣u與下三角矩陣L的逆矩陣l相乘，最終得到原 始矩陣A的逆矩陣A=uL=ul,完成整個矩陣求逆的過程。對于n 階矩陣相乘的迭代形式可表示如下：nAl j=£ ulkxlkj(10)k=j1.4實例分析4 2 15例：給定一 4階矩陣A = 8 7 2 10通過LU分解求逆矩陣A。4 8 3 6_6 8 4 9 一解：算法過程為：A* = (L x U )一1 = U 一1 乂 L" = u ' l , 第一步：求LU矩陣LU-設(shè)L000001U0011U01U02U03L10L11000U11U12U13L

6、20L21L2200110U22U23-L30L31L32L33一 000U33通過(4) (7)式可逐步進行矩陣L和U中元素的計算，如下所示:(計算L的對角)L00 = L11 =L 22=L33 = 1,(U的第一行)U 00 = a00 =4,U01 = a 01=2,U02 = a 02 = 1,U 03 =a03(L的第一列)1a10L10 一 . .-8:2, L20_ a20_ 4a30一)_ 1, L30 - ._ §U 004U004U 004(U的第二行)U 11 = a11 -L10U01 = 7 -2 2 =3,U 12 = a12 一L10U02 = 2 -

7、21 =:0,U 13 = a13 -L10U03 = 10-25=0,(L的第二列)=5,(8 - 1 2)2=1.5,1 X:1 ,L21 - | | (a21L20U 01 ) -U 11:1 /.L31(a31 - L30UU 11(U的第三行)U 22 = a22 一 L20U 02U 23 = a23 一 L20U 03(L的第三列)01 )=31x3L21 U 12L21 U 13(8-1.5 2)3 -1 1 - 2 0 = 2,6-1 5-2 0=1,115 L32(a32 一 L30U 02 一 L31U 12) = 3(4 - 1 .5 1一三 0) =1.25,U 22

8、23(U的第四行)U 33 = a33L30U 03 f L31U 13 一 L32U 23 = 9 f 1.5550 -1.251 = 0.25;3經(jīng)迭代計算，最后得到L和U矩陣為:LU分解后L矩陣= 1 .0.0000082.000000 1.000000 1.000000 2.000000 1.503000 1-6666675分解后U矩陣， 4.000000 2.0300000.0000800.0000000.0008080.0000000,00000100.0300001.0000001.2500001-0000300.0000002.00R00R9.0000000.000000 0.

9、 000800 0.030000 1.000000S.030000 0.003000 1.000000 0.250000第二步：求L和U矩陣的逆u,l-1-iu = U 1,l = L 1;(1)求u矩陣的逆0001020Uli U12u =00U22_ 000-1U03421U13030U23002U33_0005u004)1cI、00un10 00.25_ 0 04)2U2日20u03I*3u23u33由式（9）可得矩陣U的逆的各元素計算如下:（腿11U。一4(2)u”U11化1 =1 ,、 (U01u11)U 00f1_1_、1 5 c L、 c(3)u22 = = 0.5, u2 =

10、(U12u22) = (0x 0.5) = 0U 22U1131、1 5c 、u°2 = -一但0心12 U°2u22)(2 0 1 0.5) = -0.125U°04,、11八，、1、 八(4)出3 =廠=4 u23 = - (U23u33 -(1 4) = - 21 1-u13 = - (U 12u23 U3u33) = - (0 (-2) 0 4) = 0U1131 ,、1 s 八八、u03但01叫3U 02u23U°3u33)(2 0 1 (- 2) 5 4)=-4.5U°04（2）求L矩陣的逆L0010001-d111000 一111

11、001001 :=、。L1002 =21001 =1101110L20L21L2201112101120121122L°L31、2L33_1.50.66671.251130131132000 .I 133由(8)式可得L矩陣的逆的各元素計算如下(1】00 = 1,l10 = - L*00 = -220 = -(L20%0L2】10)= 3,I30 = -(L30100 + L31I10 + L32I20) = -(1.5尺 1 + 號(-2) + 仆 3)1.916667 111=1,121=-L21111 = -2,131 = -(L31111L32121) = 0.83333 1

12、22=1,132 = TL32122)= -1.25 133 =1；所以得到L和U的逆矩陣為:族畔的逆關(guān)畔： 1.000000 0.000090 0.000080 0.030000 -2.000000 1.000000 0.000000 0.000000 3.000300 -2.008000 1.000000 0.000000 -916667 0.833333 -1.250000 1.6(40000U矩陣的逆矩陣； 0.250000 -0.16666? -0.125000 -4,500000 0.000000 0-333333 0.000000 0.000000 0.000000 0.0000

13、00 0.500000 -2.0B000O 0.000000 0.000000 4.000000(3)求A的逆矩陣由式(10)可計算得到矩陣A的逆，如下:A=u l0.25 -0.166667 -0.125 -4.51111000100.3333300-2100000.5-23_210-0004 -1.916667 0.83331.25 18.833334 -3.66667 5.5 -4.5-0.66667 0.33333 005.333333 -2.66667 3-2II_-7.666667 3.33333 -54由程序計算出的結(jié)果如下：L矩陣和u矩陣乘積12 _購0阿1.0000805.0

14、000008.0000007.0000002.00000010.0000004_ 0UC3000S3-0000006.0000008-WR000R原矩陣的逆矩陣：8.833334-3.6666675.5麗明0-4.500000-0.6666670.3333330.0000000.0000005.333333-2.&G66673.000000-2”-7.6666673.333333m麗麗d4.00000ft2、C語言程序設(shè)計及測試2.1算法c程序?qū)崿F(xiàn)#include<stdio.h>#include <string.h>#define N 4void main()

15、float aNN;float LNN,UNN,outNN, out1NN;float rNN,uNN;memset( a , 0 , sizeof(a);memset( L , 0 , sizeof(L);memset( U , 0 , sizeof(U);memset( r , 0 , sizeof(r);memset( u , 0 , sizeof(u);int n=N;int k,i,j;int flag=1;float s,t;/input a matrix/ printf("ninput A=");for(i=0;i<n;i+)for(j=0;j<n

16、;j+)scanf("%f",&aij);/figure the input matrix/ printf("輸入矩陣:n");for(i=0;i<n;i+)(for (j = 0; j < n; j+)(printf("%lf ", aij);printf("n");for(j=0;j<n;j+)a0j=a0j; 計算 U 矩陣的第一for(i=1;i<n;i+)ai0=ai0/a00; 計算 L 矩陣的第1列for(k=1;k<n;k+)(for(j=k;j<n;j+

17、)(s=0;for (i=0;i<k;i+) s=s+aki*aij; 累加akj=akj-s; / 計算 U矩陣的其他元素for(i=k+1;i<n;i+)(t=0;for(j=0;j<k;j+)t=t+aij*ajk; 累加 aik=(aik-t)/akk; 計算L矩陣的其他元素for(i=0;i<n;i+)for(j=0;j<n;j+)( if(i>j)( Lij=aij; Uij=0;/ 如果 i>j ,說明行大丁列，計算矩陣的下三角部分， 得出L的值，U的為0else( Uij=aij;if(i=j) Lij=1; 否則如果i<j,說明

18、行小丁列，計算矩陣的上三角部分，得出U的值，L的為0else Lij=0;if(U11*U22*U33*U44=0)(flag=0;printf("n逆矩陣不存在");if(flag=1)(/ 求 L 和 U 矩陣的逆for (i=0;i<n;i+) /* 求矩陣 U 的逆 */ uii=1/Uii;/ 對角元素的值，直接取倒數(shù)for (k=i-1;k>=0;k-)s=0;for (j=k+1;j<=i;j+)s=s+Ukj*uji;uki=-s/Ukk;/ 迭代計算，按列倒序依次得到每一個值，for (i=0;i<n;i+) / 求矩陣 L的逆 r

19、ii=1; /對角元素的值，直接取倒 數(shù)，這里為1for (k=i+1;k<n;k+)for (j=i;j<=k-1;j+)rki=rki-Lkj*rji; 迭代計算，按列順序依次得到每一個值/繪制矩陣LU分解后的L和U 矩陣/printf("nLU 分解后 L矩陣:")；for(i=0;i<n;i+) printf("n");for(j=0;j<n;j+)printf(" %lf",Lij);printf("nLU 分解后 U矩陣:");for(i=0;i<n;i+) printf(

20、"n");for(j=0;j<n;j+)printf(" %lf",Uij);printf("n");/繪制L和U矩陣的逆矩陣 printf("nL矩陣的逆矩陣：")； for(i=0;i<n;i+) printf("n");for(j=0;j<n;j+)printf(" %lf",rij);printf("nU矩陣的逆矩陣:");for(i=0;i<n;i+) printf("n");for(j=0;j<

21、n;j+)printf(" %lf",uij);printf("n");驗證將L和U相乘，得到原矩陣printf("nL矩陣和U矩陣乘積n");for(i=0;i<n;i+)for(j=0;j<n;j+)outij=0;for(i=0;i<n;i+)for(j=0;j<n;j+)for(k=0;k<n;k+)outij+=Lik*Ukj; 2.2數(shù)據(jù)測試(1)非滿秩矩陣for(i=0;i<n;i+)for(j=0;j<n;j+)printf("%lft",outij);pr

22、intf("rn");/將r和u相乘，得到逆矩陣printf("n原矩陣的逆矩陣:n");for(i=0;i<n;i+)for(j=0;j<n;j+)out1ij=0;for(i=0;i<n;i+)for(j=0;j<n;j+)for(k=0;k<n;k+)out1ij+=uik*rkj; for(i=0;i<n;i+)for(j=0;j<n;j+)printf("%lft",out1ij);printf("rn");1、整數(shù)矩陣input A =1 3 3 4 12 3

23、4 3 6 9 2 9 7 3 6 輸入矩陣 1.000000 2,000000 1.000B98 2-000060 6.0000m9.0HB0H87.0H00H03.0000003.0000089.0000003用明麗日4.0000004.0000002-000000逆再降不存在Pregan5/ key to continue2、小數(shù)矩陣input1.2351,3752.2322.125A=1.235 2,365 4.34 5.231 2-365 4.34 5-231 6.75 5-784 0-235 4.321 2.563輸入矩陣:1.2358001.235001.3750002,2320

24、002.36500B2.365000 6.750000 4,3210004.3430004.3400005.7840002.5630005.2310005.2310000-2350002.125000逆矩陣不存在 essany key to continue(2)滿秩矩陣1整數(shù)矩陣的測試input 11=4 2 1 5B 7 2 184 8 3 &6 8 4 9拋陣：4.0000 2.Neese 1.000000 s.eetmeB.NBeeB 7.0eeeeo 2.mm luemeLaeeooe e.meo 3.000m 6.魘明明B.eera a.meee 4 .明0000 9 廊測

25、訕分解后l矩陣：1.000000 0.080006 040H0M 0.00M00 2.eeW0 1.090000 9.00H00H 0.000H0B i.eeeeee 2 皿做 i,000000 e.eeaeeo 1.500B00 1.666667 1.250008 1.00M00LU分解后虞電4.000000 2.000000 1.000000 5.006000 e.eeMee 3.000800 0.000000 虬eewe。 0.000000 土嬲g 2.0M0W 1.000000 Gk麗 138013 土BBfl瞰I 0,000000 0.250W0L矩陣的逆矩陣：1,期000 0.00

26、0000 0.000B00 0.O00Q00 -2.006000 1.906900 0.000000 0.000009 3.BB0M) -2.00000& 1.000000 0.000000 -1J16667 0,833333 -1.2500B0 1.OMN0矩陣的逆矩陣：O.250OM -0466667 -0.125M -4.5BM3B B.800000 0.333333 0.000800 0.800000 0.000000 a.000B0 0.&00000 -2.00B00B 。,網(wǎng)前明 e.roeeee e.eQ000Q 4.WMea4.咖觥LSMMLflNBNB.MM7

27、+0000ffi2.NMVM.NMB4.BN9NU顧圈3.HMBLMM皿袖觥8JMHN4伽迦由9.瞰撅L(fēng)833334-誦郵剝 BWBi-4.580000也蹣6?L333333turnLbnm5.333333-2.6t&6673.NBMII-2.K0fffiQ7,夠?1333333-JMN2小數(shù)矩陣的測試input A=1.234 2.4S 3.23 4.176 2.348 9.01 2.09® 3.14 3.18? 3.453 2-354 1.122 2-34 2-135 4*32 0-97輸入矩陣；1.234000 2.4S000S 3.230008 4.176000 2.348000 9.010000 2-098060 3.140000 3.187R06 3.4530D0 2.354003 1.122R00 2.343000 2.135000 4.320000 0,9