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

题目要求
给定一个多维矩阵，实现该矩阵的求逆运算。
1、理论分析
矩阵的一种有效而广泛应用的分解方法是矩阵的LU三角分解，将一个n阶矩阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘积。所以首先对矩阵进行三角分解，这里采用Doolittle分解，即分解为一个下三角矩阵（对角元素为1），和一个上三角矩阵的乘积。再进行相应的处理。
所以，矩阵求逆的算法流程可表述如下：

图1 矩阵求逆流程图
1)进行LU分解；
2)对分解后的L阵（下三角矩阵）和U阵（上三角矩阵）进行求逆;；
3)L阵的逆矩阵和U阵的逆矩阵相乘，即可求得原来矩阵的逆。即：
（1）
1.1矩阵的LU 分解
若n阶方阵 的各阶顺序主子式不等于零，即：
（2）
则A的LU分解存在且唯一。
（3）
由矩阵的乘法原理, 可推导出LU分解的迭代算法
（4）
（5）
（6）
（7）
矩阵的LU分解是一个循环迭代的过程, U矩阵是从第1行迭代到第n行, 而L矩阵则是从第1列迭代到第n列, 且U矩阵先于L矩阵一个节拍。
1.2 L矩阵和U矩阵求逆
首先假设下三角矩阵L的逆矩阵为，不失一般性，考虑4阶的情况，利用，有：
（1） ，，；
（2）
（3）
（4）。
从而求得下三角矩阵L的逆矩阵R式如下：
， （8）
上三角矩阵U的逆矩阵可以由下式得到：。
， （9）
矩阵求逆是一个迭代的过程，依次循环, 迭代次, 求出整个逆矩阵。其中U矩阵的循环迭代时按行顺序，列倒序进行，L矩阵的循环迭代按列顺序，行顺序进行，直到计算出整个矩阵的所有结果为止。
1.3 矩阵相乘
上三角矩阵U的逆矩阵u与下三角矩阵L的逆矩阵相乘, 最终得到原始矩阵A的逆矩阵, 完成整个矩阵求逆的过程。对于n阶矩阵相乘的迭代形式可表示如下：
（10）
1.4 实例分析
例：给定一4阶矩阵，通过LU分解求逆矩阵。
解：算法过程为：，
第一步：求LU矩阵
设，通过（4）~（7）式可逐步进行矩阵L和U中元素的计算，如下所示：

经迭代计算，最后得到L和U矩阵为：

第二步：求L和U矩阵的逆,

（1）求U矩阵的逆

由式（9）可得矩阵U的逆的各元素计算如下：

（2）求L矩阵的逆
由（8）式可得L矩阵的逆的各元素计算如下

所以得到L和U的逆矩阵为：

（3）求A的逆矩阵
由式（10）可计算得到矩阵A的逆，如下：
由程序计算出的结果如下：

2、C语言程序设计及测试
2.1 算法c程序实现

13

#include<stdio.h>
#include <string.h>
#define N 4
void main()
{ float a[N][N];
float L[N][N],U[N][N],out[N][N], out1[N][N];
float r[N][N],u[N][N];
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;j++)
scanf("%f",&a[i][j]);
//////////////////figure the input matrix//////////////////////////
printf("输入矩阵:\n");
for(i=0;i<n;i++)
{
for (j = 0; j < n; j++)
{
printf("%lf ", a[i][j]);
}
printf("\n");
}
for(j=0;j<n;j++)
a[0][j]=a[0][j]; //计算U矩阵的第一行

for(i=1;i<n;i++)
a[i][0]=a[i][0]/a[0][0]; //计算L矩阵的第1列
for(k=1;k<n;k++)
{
for(j=k;j<n;j++)
{
s=0;
for (i=0;i<k;i++)
s=s+a[k][i]*a[i][j]; //累加
a[k][j]=a[k][j]-s; //计算U矩阵的其他元素
}
for(i=k+1;i<n;i++)
{
t=0;
for(j=0;j<k;j++)
t=t+a[i][j]*a[j][k]; //累加
a[i][k]=(a[i][k]-t)/a[k][k]; //计算L矩阵的其他元素
}
}
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{ if(i>j)
{ L[i][j]=a[i][j]; U[i][j]=0;}//如果i>j，说明行大于列，计算矩阵的下三角部分，得出L的值，U的//为0
else
{ U[i][j]=a[i][j];
if(i==j) L[i][j]=1; //否则如果i<j，说明行小于列，计算矩阵的上三角部分，得出U的//值，L的为0
else L[i][j]=0;
}
} if(U[1][1]*U[2][2]*U[3][3]*U[4][4]==0){
flag=0;
printf("\n逆矩阵不存在");}
if(flag==1){
/////////////////////求L和U矩阵的逆
for (i=0;i<n;i++) /*求矩阵U的逆 */
{u[i][i]=1/U[i][i];//对角元素的值，直接取倒数
for (k=i-1;k>=0;k--)
{s=0;
for (j=k+1;j<=i;j++)
s=s+U[k][j]*u[j][i];
u[k][i]=-s/U[k][k];//迭代计算，按列倒序依次得到每一个值，
}
}
for (i=0;i<n;i++) //求矩阵L的逆
{r[i][i]=1; //对角元素的值，直接取倒数，这里为1
for (k=i+1;k<n;k++)
{for (j=i;j<=k-1;j++)
r[k][i]=r[k][i]-L[k][j]*r[j][i]; //迭代计算，按列顺序依次得到每一个值
}
}
/////////////////绘制矩阵LU分解后的L和U矩阵///////////////////////
printf("\nLU分解后L矩阵:");
for(i=0;i<n;i++)
{ printf("\n");
for(j=0;j<n;j++)
printf(" %lf",L[i][j]);
}
printf("\nLU分解后U矩阵:");
for(i=0;i<n;i++)
{ printf("\n");
for(j=0;j<n;j++)
printf(" %lf",U[i][j]);
}
printf("\n");
////////绘制L和U矩阵的逆矩阵
printf("\nL矩阵的逆矩阵:");
for(i=0;i<n;i++)
{ printf("\n");
for(j=0;j<n;j++)
printf(" %lf",r[i][j]);
}
printf("\nU矩阵的逆矩阵:");
for(i=0;i<n;i++)
{ printf("\n");
for(j=0;j<n;j++)
printf(" %lf",u[i][j]);
}
printf("\n");
//验证将L和U相乘，得到原矩阵
printf("\nL矩阵和U矩阵乘积\n");
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{out[i][j]=0;}
}
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
for(k=0;k<n;k++)
{out[i][j]+=L[i][k]*U[k][j];}
}
}
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
printf("%lf\t",out[i][j]);
}
printf("\r\n");
}
//////////将r和u相乘，得到逆矩阵
printf("\n原矩阵的逆矩阵:\n");
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{out1[i][j]=0;}
}
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
for(k=0;k<n;k++)
{out1[i][j]+=u[i][k]*r[k][j];}
}
}
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
printf("%lf\t",out1[i][j]);
}
printf("\r\n");
}
}
}

2.2 数据测试
（1）非满秩矩阵
1>、整数矩阵

2>、小数矩阵

（2）满秩矩阵
1> 整数矩阵的测试

2> 小数矩阵的测试