矩阵LU分解求逆详细分析与C语言实现.docx
- 文档编号:1491922
- 上传时间:2022-10-22
- 格式:DOCX
- 页数:15
- 大小:134.04KB
矩阵LU分解求逆详细分析与C语言实现.docx
《矩阵LU分解求逆详细分析与C语言实现.docx》由会员分享,可在线阅读,更多相关《矩阵LU分解求逆详细分析与C语言实现.docx(15页珍藏版)》请在冰豆网上搜索。
矩阵LU分解求逆详细分析与C语言实现
题目要求
给定一个多维矩阵,实现该矩阵的求逆运算。
1、理论分析
矩阵的一种有效而广泛应用的分解方法是矩阵的LU三角分解,
将一个n阶矩阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘
积。
所以首先对矩阵进行三角分解,这里采用Doolittle分解,即分解为一个下三角矩阵(对角元素为1),和一个上三角矩阵的乘积。
再进行相应的处理。
所以,矩阵求逆的算法流程可表述如下:
图1矩阵求逆流程图
1)进行LU分解;
2)对分解后的L阵(下三角矩阵)和U阵(上三角矩阵)进行求逆;;
3)L阵的逆矩阵和
U阵的逆矩阵相乘,即可求得原来矩阵的逆。
即:
A1(LU)
1U1L1
1.1矩阵的LU
分解
若n阶方阵A
C:
n的各阶顺序主子式不等于零,即:
a11
a12
a1k
k
a21
a22
a2k
0,(k
1,2,
n),
ak1
ak2
akk
则A的LU分解A
L
U存在且唯一。
a11
L
a1r
L
a1n
M
M
M
A
ar1
L
arr
L
arn
M
M
M
an1
L
anr
L
ann
1
U11
L
U1r
L
U1n
M
O
O
M
Lr1
L
1
Urr
L
Urn
M
M
O
O
M
Ln1
K
Lnr
L
1
Unn
(3)
LU
由矩阵的乘法原理,可推导出
LU分解的迭代算法
Uoj
aoj,(j0,1,2,L,n
1),
li0
Urj
a
宀(i0,1,2,L,n1),
Uoo
r1
arjlikUkj
Irr
arj
r1
1rkukj,
k1
(r
0,1,2,L,n1;jr,L,n1),
r1
airlikUkr
lir
(r
k1
Urr
0,1,2,L,n1;ir1,L,n1)
矩阵的LU分解是一个循环迭代的过程,U矩阵是从第1行迭代到第n
行,而L矩阵则是从第1列迭代到第n列,且U矩阵先于L矩阵一个
节拍。
1.2L矩阵和U矩阵求逆
首先假设下三角矩阵L的逆矩阵为I,不失一般性,考虑4阶的情况,
「30
r00(r31L-iQr32L2Or33L30)L。
从而求得下三角矩阵L的逆矩阵R式如下:
lii1i
j,i
rji
ri(rjklki),iki1.
0,i
上三角矩阵U的逆矩阵可以由下式得到:
。
矩阵求逆是一个迭代的过程,依次循环,迭代n1次,求出整个
逆矩阵。
其中U矩阵的循环迭代时按行顺序,列倒序进行,L矩阵的
循环迭代按列顺序,行顺序进行,直到计算出整个矩阵的所有结果为
止。
1.3矩阵相乘
上三角矩阵U的逆矩阵u与下三角矩阵L的逆矩阵I相乘,最终得到原
始矩阵A的逆矩阵A1U1L1ul,完成整个矩阵求逆的过程。
对于n
阶矩阵相乘的迭代形式可表示如下:
第一步:
求LU矩阵
*-00
22200
01
202
203
L10
L11
0
0
0
U11
U12
U13
L20
L21
L22
0
0
0
U22
U23
L30
L31
L32
L33
0
0
0
U33
通过(4)
-(7)式可逐步进行矩阵L和U中元素的计算,如下所示:
(计算L的对角)
L00L11
(U的第一行
U00a00
(L的第一列
a10
L22L331,
)
4,U01
)
8
4
)
L10U
a01
2,U02
a021,U03a035,
2,
L20
a20
L10U
U00(U的第二行
U11a11
U12a12
U13a13
(L的第二列
1
-(a21L20U01)
U11
1
U(a31L30U01)
U11
(U的第三行)
U22a22L20U02
U23*23L20U03
(L的第三列)
1
L32(a32L30U02L31U12)
U22
(U的第四行)
1丄30
a30
U00
1.5,
01
L10U02
L10U03
)
7
2
10
3,
0,
0,
L21
L31
1
3
1
3
(8
(8
2)
1.5
2)
L21U12
L21U13
U33a33L30U03L31U13L32U23
经迭代计算,最后得到L和U矩阵为:
(4
02,
01,
1.51
1.5
0)1.25,
01.2510.25;
第二步:
求L和U矩阵的逆u,l
uU1,lL1;
(1)求U矩阵的逆
U00
U01
U02
U03
4
2
1
5
u00
u01
u02
u03
0
U11
U12
U13
0
3
0
0
0
u11
u12
u13
0
0
U22
U23
0
0
2
1
0
0
u22
u23
0
0
0
U33
0
0
0
0.25
0
0
0
u33
求L矩阵的逆
1-00
0
0
01
1
0
0
01
l000
0
0
L10
L11
0
0
2
1
0
0
l10l11
0
0
I
L20
L21
L22
0
1
2
1
0
l20l21
l22
0
L30
L31
L32
L33
1.5
0.666
71.25
1
l30l31
l32
l33
⑷l331;
所以得到L和U的逆矩阵为:
求A的逆矩阵
由式(10)可计算得到矩阵A的逆,如下:
A1
由程序计算出的结果如下:
2、C语言程序设计及测试2.1算法c程序实现
#includevstdio.h>
#inelude
#defineN4
voidmain()
{floata[N][N];
floatL[N][N],U[N][N],out[N][N],out1[N][N];
floatr[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));
intn=N;
intk,i,j;
intflag=1;
floats,t;
////////////////////inputamatrix////printf("\ninputA=");
for(i=0;i for(j=0;j //////////////////figuretheinput matrix////////////////////////// printf("输入矩阵: \n"); for(i=0;i { for(j=0;j { printf("%lf",a[i][j]); } printf("\n"); } for(j=0;j a[0][j]=a[0][j];〃计算U矩阵的第一 行 for(i=1;ivn;i++) a[i][0]=a[i][0]/a[0][0];//计算L矩 阵的第1列 for(k=1;k { for(j=k;j { s=0; for(i=0;i 〃累加 //计算 s=s+a[k][i]*a[i][j];〃累加 a[k][j]=a[k][j]-s;//计算U矩阵的其他元素 } for(i=k+1;ivn;i++) { t=O; for(j=O;jvk;j++) t=t+a[i][j]*a[j][k];a[i][k]=(a[i][k]-t)/a[k][k]; L矩阵的其他元素 } } for(i=O;ivn;i++) for(j=O;jvn;j++) {if(i>j) {L[i][j]=a[i][j];U[i][j]=O;}//如果i>j,说明行大于列,计算矩阵的下三角部分,得出L的值,U的//为O else {U[i][j]=a[i][j];if(i==j)L[i][j]=1;//否则如果 ivj,说明行小于列,计算矩阵的上角部分,得出U的〃值,L的为O elseL[i][j]=O; } }if(U[1][1]*U[2][2]*U[3][3]*U[4][4]==O){ flag=O; printf("\n逆矩阵不存在");}if(flag==1){/////////////////////求L和U矩阵的逆 for(i=O;ivn;i++)/*求矩阵U的逆*/{u[i][i]=1/U[i][i];//对角元素的值,直接取倒数 for(k=i-1;k>=O;k--) {s=O; for(j=k+1;jv=i;j++)s=s+U[k][j]*u[j][i];u[k][i]=-s/U[k][k];//迭代计算,按列倒序依次得到每一个值, } } for(i=O;ivn;i++)//求矩阵L的逆{r[i][i]=1;//对角元素的值,直接取倒数,这里为1 for(k=i+1;kvn;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;ivn;i++) {printf("\n"); for(j=O;jvn;j++) printf("%lf",L[i][j]); } printf("\nLU分解后U矩阵: "); for(i=O;ivn;i++) {printf("\n"); for(j=O;jvn;j++) printf("%lf",U[i][j]); } printf("\n"); ////////绘制L和U矩阵的逆矩阵 printf("\nL矩阵的逆矩阵: "); for(i=O;i
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵 LU 分解 详细 分析 语言 实现