矩阵LU分解求逆详细分析与C语言实现.docx
题目要求给定一个多维矩阵，实现该矩阵的求逆运算。1、理论分析矩阵的一种有效而广泛应用的分解方法是矩阵的1.U三角分解，将一个n阶矩阵A分解为一个下三角矩阵1.和一个上三角矩阵U的乘积。所以首先对矩阵进行三角分解,这里采用DOolittIe分解，即分解为一个下三角矩阵（对角元素为1）,和一个上三角矩阵的乘积。再进行相应的处理。所以，矩阵求逆的算法流程可表述如下：图1矩阵求逆流程图1）进行1.U分解；2）对分解后的1.阵（下三角矩阵）和U阵（上三角矩阵）进行求逆;；3）1.阵的逆矩阵和U阵的逆矩阵相乘，即可求得原来矩阵的逆。即：A-,=（1.U）T=UT1.（1）1.1矩阵的1.U分解假设n阶方阵的各阶顺序主子式不等于零，即：a1¢/12aka2a22a2kr.a=.0,（%=1,2,），（2）akak2akk那么A的1.U分解A=1.XU存在且唯一。由矩阵的乘法原理，可推导出1.U分解的迭代算法UOj=4)/,(/=0,1,2,.,h-1),IiO=i-0,l,2,11-l),wOOr-%-EIlikUkjZ-I%=%一*%,Irrk=l(r=0,1,2,H1;7=r.,H1),r-lairEjiMkr1 k=l(7)(r=0,1,2,1,111；i=r+1,/21)矩阵的1.U分解是一个循环迭代的过程,U矩阵是从第1行迭代到第n行，而1.矩阵那么是从第1列迭代到第n歹J,且U矩阵先于1.矩阵一个节拍。1.2 1.矩阵和U矩阵求逆首先假设下三角矩阵1.的逆矩阵为I,不失一般性，考虑4阶的情况,利用=/,有:(1)100=4)1.l=耳，,22=5；，,33=4；»4o=T0O(411.lo);(3)，20=，00(，21乙10+，22心20);4Go=-a(4140+4240+公及)。从而求得下三角矩阵1.的逆矩阵R式如下:/Tjg7=1-(),<7上三角矩阵U的逆矩阵可以由下式得到：O1i'I=J一(EUjkUki),i>j%5,i<j0矩阵求逆是一个迭代的过程，依次循环，迭代-1次，求出整个逆矩阵。其中U矩阵的循环迭代时按行顺序，列倒序进行，1.矩阵的循环迭代按列顺序,行顺序进行，直到计算出整个矩阵的所有结果为止。1.3 矩阵相乘上三角矩阵U的逆矩阵U与下三角矩阵1.的逆矩阵/相乘，最终得到原始矩阵A的逆矩阵AT=UT1.T=R,完成整个矩阵求逆的过程。对于n阶矩阵相乘的迭代形式可表示如下：(10)*皿刃团伙x伙刃k=j1.4实例分析4215通过1.U分解求逆矩阵A-'例：给定一4阶矩阵A=：；?48366849解:算法过程为：A-=(1.XST=UX1.=IIX第一步：求1.U矩阵-1.oOo0。一UOOUOl002U031.XU=J1.IO1.Il0°X0UnUnU3设乙201.?Ii22。00U22U23，通过(4)_乙304243_OOOU33(7)式可逐步进行矩阵1.和U中元素的计算，如下所示:经迭代计算，最后得到1.和U矩阵为：第二步：求1.和U矩阵的逆,/(1)求U矩阵的逆由式(9)可得矩阵U的逆的各元素计算如下:(2)求1.矩阵的逆4k一.-000-1-1000"-I,00000A1002100-I，10I1100'2201210V,20220由11.32%1.50.66671.251_,3041Iyl/33_式可得1.矩阵的逆的各元素计算如下所以得到1.和U的逆矩阵为：求A的逆矩阵由式(10)可计算得到矩阵A的逆，如下:A-1=u×l0.25-0.166667-0.125-4.5'1000"00.3333300×-2100000.5-23-2100004_-1.9166670.83333-1.251.由程序8.833334-3.666675.5-4.5-0.666670.33333005.333333-2.666673-2-7.6666673.33333-54计算出的结果如下:2、C语言程序设计及测试2.1算法C程序实现#include<stdio.h>#include<string.h>#defineN4voidmain()floataNN;float1.NN,UNN,outNN,outUNN;floatrNN,uNN;memset(a,0,sizeof（八）);memset(1.,0,sizeof(1.);memset(U,0,sizeof(U);memset(r,0,sizeof(r);memset(u,0,sizeof(u);intn=N;intk,i,j;intflag=1;floats,t;/inputamatrix/printf(,'ninputA=m);for(i=0;i<n;i+)for(j=0;j<n;j+)scanf(,%f&aij);/figuretheinputmatrix/Printf("输入矩阵:n)for(i=0;i<n;i+)(for(j=0;j<n;j+)(p11ntf(11%lf11,alij);)printf(11nn);)for(j=0;j<n;j+)alj=alj;计算U矩阵的第一行for(i=l;i<n;i+)ai0=ai0a00;计算1.矩阵的第1列for(k=l;k<n;k+)(for(j=k;j<n;j+)(s=0;for(i=0;i<k;i+)s=s+aki*aig;累加ak11j=akj-s;计算U矩阵的其他元素for(i=k+l;i<n;i+)(t=0;for(j=0;j<k;j+)t=t+aij*ajk;累加aik=(aik-t)akk;计算1.矩阵的其他元素)fdr(i=0;i<n;i+)for(j=0;j<n;j+)if(i>j)皿=0;如果i>j,说明行大于列，计算矩阵的下三角局部,得出1.的值，U的为0elseUifj=aij;if(i=j)1.ij=l;否那么如果ivj，说明行小于列，计算矩阵的上三角局部，得出U的/值，1.的为0else1.ij=0;if(Ull*U22*U33*U44=0)(flag=0;Printf("n逆矩阵不存在");if(flag=l)/求1.和U矩阵的逆for(i=0;ivn;i+)/*求矩阵U的逆*/uii=lUii;对角元素的值,直接取倒数for(k=i-l;k>=0;k-)s=0;for(j=k+hj<=iy+÷)s=s+Ukj*uji;uk11i=sUkk;迭代计算，按列倒序依次得到每一个值，for(i=0;i<n;i+)/求矩阵1.的逆rii=l;对角元素的值，直接取倒数，这里为1for(k=i+l;k<n;k+)for(j=i;j<=k-l;j+)rki=rki-1.kj*rji;迭代计算，按列顺序依次得到每一个值/绘制矩阵1.U分解后的1.和U矩阵/PrindrW1.U分解后1.矩阵:")；for(i=0;i<n;i+)printf(11n");for(j=0;j<n;j+)printf(,%l,1.ij);)PrindrW1.U分解后U矩阵:")；for(i=0;i<n;i+)printf(',n11);for(j=0;j<n;j+)P1111tf(11%lf,Uij);)printf(,n11);/绘制1.和U矩阵的逆矩阵Printfn1.矩阵的逆矩阵:")；for(i=0;i<n;i+)printf(',n');for(j=0;j<n;j+)printf(11%l,rij);)Printf(11U矩阵的逆矩阵巧;fbr(i=0;i<n;i+)printf(,n*');for(j=0;j<n;j+)P11ntf(,%lf,uij);Iprintf(11nu);验证将1.和U相乘，得到原矩阵Printf(n1.矩阵和U矩阵乘积n");for(i=0;i<n;i+)forG=0;j<n;j+)outij=0;)for(i=0;i<n;i+)(for(j=0;j<n;j+)(for(k=0;k<n;k+)outifj÷=1.ik*Ukj;)for(i=0;i<n;i+)(for(j=0;j<n;j+)(printf("%lt",outiirj);Prindrrn")；)/将和U相乘，得到逆矩阵Printlrn原矩阵的逆矩阵:n)for(i=0;i<n;i+)2.2数据测试非满秩矩阵1>、整数矩阵2>、小数矩阵(2)满秩矩阵1>整数矩阵的测试2>小数矩阵的测试for(j=0;j<n;j+)outlij=0;)for(i=0;i<n;i+)(for(j=0;j<n;j+)(for(k=0;k<n;k+)outlij+=uik*rkj;)for(i=0;i<n;i+)(for(j=0;j<n;j+)(printf(11%lt11,outlirj);)printf(11rn11);)
