矩阵LU分解求逆详细分析与C语言实现

1、题目要求给定一个多维矩阵，实现该矩阵的求逆运算。1、理论分析矩阵的一种有效而广泛应用的分解方法是矩阵的 LU三角分解， 将一个n阶矩阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘 积。所以首先对矩阵进行三角分解，这里采用 Doolittle分解，即分解 为一个下三角矩阵（对角元素为1）,和一个上三角矩阵的乘积。再进 行相应的处理。所以，矩阵求逆的算法流程可表述如下：'输输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 -由矩阵的乘法原理,可推导出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分解是一个循环迭代的过程，U矩阵是从第1行迭代到第n 行,而L矩阵则是从第1列迭代到第n列，且U矩阵先于L矩阵一个 节拍。1.2 L矩阵和U矩阵求逆首先假设下三角矩阵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矩阵求逆是一个迭代的过程，依次循环，迭代n1次，求出整个逆矩阵。其中U矩阵的循环迭代时按行顺序，列倒序进行， L矩阵的 循环迭代按列顺序，行顺序进行，直到计算

5、出整个矩阵的所有结果为 止。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-设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经迭代计算，最后得到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的逆矩阵为:族畔的逆关畔： 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由程序计算出的结果如下：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语言程序设计及测试2.1算法c程序实现#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;/ 对角元素的值，直接取倒数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; /对角元素的值，直接取倒 数，这里为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数据测试(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、整数矩阵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、小数矩阵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整数矩阵的测试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 。,网前明 e.roeeee e.eQ000Q 4.WMea4.咖觥LSMMLflNBNB.MM7

27、+0000ffi2.NMVM.NMB4.BN9NU顾圈3.HMBLMM皿袖觥8JMHN4伽迦由9.瞰撅L833334-诵邮剥 BWBi-4.580000也蹒6?L333333turnLbnm5.333333-2.6t&6673.NBMII-2.K0fffiQ7,够?1333333-JMN2小数矩阵的测试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