矩陣LU分解求逆詳細(xì)分析與C語(yǔ)言實(shí)現(xiàn)

1、題目要求給定一個(gè)多維矩陣，實(shí)現(xiàn)該矩陣的求逆運(yùn)算。1、理論分析矩陣的一種有效而廣泛應(yīng)用的分解方法是矩陣的LU三角分解，將一個(gè)n階矩陣A分解為一個(gè)下三角矩陣L和一個(gè)上三角矩陣U的乘積。所以首先對(duì)矩陣進(jìn)行三角分解，這里采用Doolittle分解，即分解為一個(gè)下三角矩陣（對(duì)角元素為1），和一個(gè)上三角矩陣的乘積。再進(jìn)行相應(yīng)的處理。所以，矩陣求逆的算法流程可表述如下：圖1 矩陣求逆流程圖1)進(jìn)行LU分解；2)對(duì)分解后的L陣（下三角矩陣）和U陣（上三角矩陣）進(jìn)行求逆;；3)L陣的逆矩陣和U陣的逆矩陣相乘，即可求得原來(lái)矩陣的逆。即： （1）1.1矩陣的LU 分解若n階方陣 的各階順序主子式不等于零，即： （2

2、）則A的LU分解存在且唯一。 （3）由矩陣的乘法原理, 可推導(dǎo)出LU分解的迭代算法（4） （5） （6） （7）矩陣的LU分解是一個(gè)循環(huán)迭代的過(guò)程, U矩陣是從第1行迭代到第n行, 而L矩陣則是從第1列迭代到第n列, 且U矩陣先于L矩陣一個(gè)節(jié)拍。1.2 L矩陣和U矩陣求逆首先假設(shè)下三角矩陣L的逆矩陣為，不失一般性，考慮4階的情況，利用，有：（1） ，；（2）（3）（4）。從而求得下三角矩陣L的逆矩陣R式如下：， （8）上三角矩陣U的逆矩陣可以由下式得到：。， （9）矩陣求逆是一個(gè)迭代的過(guò)程，依次循環(huán), 迭代次, 求出整個(gè)逆矩陣。其中U矩陣的循環(huán)迭代時(shí)按行順序，列倒序進(jìn)行，L矩陣的循環(huán)迭代按列順

3、序，行順序進(jìn)行，直到計(jì)算出整個(gè)矩陣的所有結(jié)果為止。1.3 矩陣相乘上三角矩陣U的逆矩陣u與下三角矩陣L的逆矩陣相乘, 最終得到原始矩陣A的逆矩陣, 完成整個(gè)矩陣求逆的過(guò)程。對(duì)于n階矩陣相乘的迭代形式可表示如下： （10）1.4 實(shí)例分析例：給定一4階矩陣，通過(guò)LU分解求逆矩陣。解：算法過(guò)程為：，第一步：求LU矩陣設(shè)，通過(guò)（4）（7）式可逐步進(jìn)行矩陣L和U中元素的計(jì)算，如下所示：經(jīng)迭代計(jì)算，最后得到L和U矩陣為：第二步：求L和U矩陣的逆,（1）求U矩陣的逆由式（9）可得矩陣U的逆的各元素計(jì)算如下：（2）求L矩陣的逆由（8）式可得L矩陣的逆的各元素計(jì)算如下所以得到L和U的逆矩陣為：（3）求A的逆矩

4、陣由式（10）可計(jì)算得到矩陣A的逆，如下：由程序計(jì)算出的結(jié)果如下：2、C語(yǔ)言程序設(shè)計(jì)及測(cè)試2.1 算法c程序?qū)崿F(xiàn)13#include#include #define N 4void main() 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

5、 k,i,j; int flag=1; float s,t; /input a matrix/ printf(ninput A=); for(i=0;in;i+) for(j=0;jn;j+) scanf(%f,&aij);/figure the input matrix/printf(輸入矩陣:n);for(i=0;in;i+)for (j = 0; j n; j+)printf(%lf , aij);printf(n);for(j=0;jn;j+) a0j=a0j; /計(jì)算U矩陣的第一行 for(i=1;in;i+) ai0=ai0/a00; /計(jì)算L矩陣的第1列 for(k=1;kn;k

6、+) for(j=k;jn;j+) s=0; for (i=0;ik;i+) s=s+aki*aij; /累加 akj=akj-s; /計(jì)算U矩陣的其他元素 for(i=k+1;in;i+) t=0; for(j=0;jk;j+) t=t+aij*ajk; /累加 aik=(aik-t)/akk; /計(jì)算L矩陣的其他元素 for(i=0;in;i+) for(j=0;jj) Lij=aij; Uij=0;/如果ij，說(shuō)明行大于列，計(jì)算矩陣的下三角部分，得出L的值，U的/為0 else Uij=aij; if(i=j) Lij=1; /否則如果ij，說(shuō)明行小于列，計(jì)算矩陣的上三角部分，得出U的/

7、值，L的為0 else Lij=0; if(U11*U22*U33*U44=0) flag=0; printf(n逆矩陣不存在); if(flag=1) /求L和U矩陣的逆 for (i=0;i=0;k-)s=0;for (j=k+1;j=i;j+)s=s+Ukj*uji;uki=-s/Ukk;/迭代計(jì)算，按列倒序依次得到每一個(gè)值，for (i=0;in;i+) /求矩陣L的逆 rii=1; /對(duì)角元素的值，直接取倒數(shù)，這里為1for (k=i+1;kn;k+)for (j=i;j=k-1;j+)rki=rki-Lkj*rji; /迭代計(jì)算，按列順序依次得到每一個(gè)值/繪制矩陣LU分解后的L和U

8、矩陣/ printf(nLU分解后L矩陣:); for(i=0;in;i+) printf(n); for(j=0;jn;j+) printf( %lf,Lij); printf(nLU分解后U矩陣:); for(i=0;in;i+) printf(n); for(j=0;jn;j+) printf( %lf,Uij); printf(n); /繪制L和U矩陣的逆矩陣 printf(nL矩陣的逆矩陣:); for(i=0;in;i+) printf(n); for(j=0;jn;j+) printf( %lf,rij); printf(nU矩陣的逆矩陣:); for(i=0;in;i+) pr

9、intf(n); for(j=0;jn;j+) printf( %lf,uij); printf(n); /驗(yàn)證將L和U相乘，得到原矩陣 printf(nL矩陣和U矩陣乘積n); for(i=0;in;i+) for(j=0;jn;j+)outij=0; for(i=0;in;i+) for(j=0;jn;j+) for(k=0;kn;k+) outij+=Lik*Ukj; for(i=0;in;i+) for(j=0;jn;j+) printf(%lft,outij); printf(rn); /將r和u相乘，得到逆矩陣 printf(n原矩陣的逆矩陣:n); for(i=0;in;i+) for(j