LU和QR分解法解线性方程组.docx
- 文档编号:22785073
- 上传时间:2023-04-27
- 格式:DOCX
- 页数:15
- 大小:46.71KB
LU和QR分解法解线性方程组.docx
《LU和QR分解法解线性方程组.docx》由会员分享,可在线阅读,更多相关《LU和QR分解法解线性方程组.docx(15页珍藏版)》请在冰豆网上搜索。
LU和QR分解法解线性方程组
LU和QR法解线性方程组
一、问题描述
求解方程组
==
要求:
1、编写用三角〔LU〕分解法求解线性方程组;
2、编写用正交三角〔QR〕分解法求解线性方程组.
二、问题分析
求解线性方程组Ax=b,其实质就是把它的系数矩阵A通过各种变换成一个下三角或上三角矩阵,从而简化方程组的求解.因此,在求解线性方程组的过程中,把系数矩阵A变换成上三角或下三角矩阵显得尤为重要,然而矩阵A的变换通常有两种分解方法:
LU分解法和QR分解法.
1、LU分解法:
将A分解为一个下三角矩阵L和一个上三角矩阵U,即:
A=LU,
其中L=
U=
2、QR分解法:
将A分解为一个正交矩阵Q和一个上三角矩阵R,即:
A=QR
三、实验原理
1、LU分解法
解Ax=b的问题就等价于要求解两个三角形方程组:
⑴Ly=b,求y;
⑵Ux=y,求x.
设A为非奇异矩阵,且有分解式A=LU,L为单位下三角阵,U为上三角阵.
L,U的元素可以有n步直接计算定出.用直接三角分解法解Ax=b〔要求A的所有顺序主子式都不为零〕的计算公式:
①
i=2,3,…,n.
计算U的第r行,L的第r列元素:
②
i=r,r+1,…,n;
③
i=r+1,…,n,且r≠n.
求解Ly=b,Ux=y的计算公式;
④
⑤
2、QR分解法
四、实验步骤
1、LU分解法
1>将矩阵A保存进计算机中,再定义2个空矩阵L,U以便保存求出的三角矩阵的值.
利用公式①,②,③将矩阵A分解为LU,L为单位下三角阵,U为上三角阵.
2>可知计算方法有三层循环.
先通过公式①计算出U矩阵的第一行元素
和L矩阵的第一列元素
.
再根据公式②和③,和上次的出的值,求出矩阵其余的元素,每次都要三次循环,求下一个元素需要上一个结果.
3>先由公式④,Ly=b
求出y,因为L为下三角矩阵,所以由第一行开始求y.
4>再由公式⑤,Ux=y
求出x,因为U为上三角矩阵,所以由最后一行开始求x.
2、QR分解法
五、程序流程图
1、LU分解法
2、QR分解法
六、实验结果
1、LU分解法
2、QR分解法
七、实验总结
为了求解线性方程组,我们通常需要一定的解法.其中一种解法就是通过矩阵的三角分解来实现的,属于求解线性方程组的直接法.在不考虑舍入误差下,直接法可以用有限的运算得到精确解,因此主要适用于求解中小型稠密的线性方程组.
1、三角分解法
三角分解法是将A矩阵分解成一个上三角形矩阵U和一个下三角形矩阵L,这样的分解法又称为LU分解法.它的用途主要在简化一个大矩阵的行列式值的计算过程,求反矩阵和求解联立方程组.不过要注意这种分解法所得到的上下三角形矩阵并非唯一,还可找到数个不同的一对上下三角形矩阵,此两三角形矩阵相乘也会得到原矩阵.
2、QR分解法
QR分解法是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为QR分解法.
在编写这两个程序过程中,起初遇到不少麻烦!
虽然课上老师反复重复着:
"算法不难的,It'ssoeasy!
"但是当自己实际操作时,感觉并不是那么容易.毕竟是要把实际的数学问题转化为计算机能够识别的编程算法,所以在编写程序之前我们仔细认真的把所求解的问题逐一进行详细的分析,最终转化为程序段.每当遇到问题时,大家或许有些郁闷,但最终还是静下心来反复仔细的琢磨,一一排除了错误,最终完成了本次实验.
回头一想原来编个程序其实也没有想象的那么复杂,只要思路清晰,逐步分析,就可以慢慢搞定了.
通过这次实验,让我们认知到团队的作用力度是不容忽视的,以后不管干任何时都要注重团队合作,遇到不懂得不明白的大家一起讨论,越讨论越清晰,愈接近最优答案.这样不管干什么都能事半功倍.庆幸自己有这么个团队,也明白大家一起劳动的果实最珍贵.
附源代码:
LR分解法:
#include
voidinput_A<>;//输入矩阵A
voidinput_b<>;//输入矩阵b
voidoutput_x
voidmain<>
{
inti,k,r;
floatA[4][4]={{4,2,1,5},{8,7,2,10},{4,8,3,6},{12,6,11,20}},b[4]={-2,-7,-7,-3},x[4],u[4][4],l[4][4],y[4];//给定的方程组
//input_A<>;
//input_b<>;
//显示矩阵A、b
printf<"矩阵A[4][4]:
\n">;
for
{
for
printf<"%-10f",A[i][k]>;
printf<"\n">;
}
for
u[0][i]=A[0][i];
for
{
l[i][0]=A[i][0]/u[0][0];
l[i][i]=1;
}
for
{
for
{
floatt=0;
for
t=t+l[r][k]*u[k][i];
u[r][i]=A[r][i]-t;
}
for
{
floatt1=0;
for
t1+=l[i][k]*u[k][r];
l[i][r]=/u[r][r];
}
}
y[0]=b[0];
for
{
floatt=0;
for
t+=l[i][k]*y[k];
y[i]=b[i]-t;
}
for=0;i-->
{
floatt=0;
for
t+=u[i][k]*x[k];
x[i]=
}
printf<"*****************************\n">;
output_x
}
//输入矩阵A
voidinput_A<>
{
floatA[4][4];
inti,j;
printf<"inputmatrixA[4][4]:
\n">;
for
for
scanf<"%f",&A[i][j]>;
}
//输入矩阵b
voidinput_b<>
{
floatb[4];
inti;
printf<"inputmatrixb[4]:
\n">;
for
scanf<"%f",&b[i]>;
}
//输出方程的根x
voidoutput_x
{
inti;
printf<"方程组的根为:
\n">;
for
printf<"x[%d]=%-13f",i+1,x[i]>;
printf<"\n">;
}
QR分解法:
#include
#include
#defineN4
voidmatrix_time
voidmatrix_vec
doublevec_value
voidvec_time
voidhouseholder
voidmatrix_turn
voidsolution
voidQR
voidmain<>
{
doubleA[4][4]={{4,2,1,5},{8,7,2,10},{4,8,3,6},{12,6,11,20}};
doubleb[4]={-2,-7,-7,-3};
inti,j;
intn=4;
printf<"待求解的方程组系数矩阵A为:
\n">;
for//显示矩阵A
{
for
printf<"%-10f",A[i][j]>;
printf<"\n">;
}
QR;
}
voidmatrix_time
{
inti,j,k;
for
for
{
C[i][j]=0;
for
C[i][j]=C[i][j]+A[i][k]*B[k][j];
}
}
voidmatrix_vec
{
inti,j;
for
{
C[i]=0;
for
C[i]=C[i]+A[i][j]*B[j];
}
}
doublevec_value
{
doubleSum=0;
inti;
for
Sum=Sum+A[i]*A[i];
Sum=sqrt
returnSum;
}
voidvec_time
{
inti,j;
for
for
H[i][j]=a[i]*a[j];
}
voidhouseholder
{
inti,j;
doublefirst;//存放首元素
doublevec_v;//存放向量的模
doublea1[N]={0};
for
for
{
if
H[i][j]=1;
else
H[i][j]=0;
}
first=a[m];//取矩阵A转置的第m行〔取矩阵A的第m列数〕
vec_v=vec_value<&a[m],n-m>;//计算m列元素所构成向量的模
a[m]=a[m]-vec_v;//w的分子部分
vec_v=vec_value<&a[m],n-m>;//w的分母部分
vec_v=vec_v*vec_v;
for
for
H[i][j]+=-a[i]*a[j]*2/vec_v;
}
voidmatrix_turn
{
doublea[N][N]={0};
inti,j;
for
for
a[i][j]=A[i][j];
for
for
A[i][j]=a[j][i];
}
voidsolution
{
inti,j;
doublex[N]={0};
doublesum=0;
for=0;i-->
{
for
sum+=A[i][j]*x[j];
x[i]=/A[i][i];
sum=0;
}
printf<"矩阵的QR分解求解结果为:
\n">;
for
printf<"X[%d]=%-10f\n",i+1,x[i]>;
}
voidQR
{
inti,j,k,t;
doubleH1[N][N]={0},H2[N][N]={0},H3[N][N]={0};//H1,H2存放相邻两次的Householder矩阵,H3为中间变量矩阵
doubletemp[N][N]={0};
doubleQ[N][N]={0},R[N][N]={0},Q1[N][N]={0};
doubleb1[N]={0};
for//将矩阵A临时存放在temp中
for
temp[i][j]=A[i][j];
for
{
H2[i][i]=1;//单位阵
}
matrix_turn
for//Householder的一次变换
{
householder
matrix_time
;//矩阵H1和A相乘放在Q中
for
for
A[k][t]=Q[k][t];
for
for
temp[k][t]=A[k][t];
matrix_turn
matrix_time
;
for
for
H2[k][t]=H3[k][t];//H2是一次变换与A相乘后的矩阵〔H1*A〕
}
for
for
R[k][t]=A[k][t];
printf<"*********************************\n">;
printf<"系数矩阵A进行QR分解后的R矩阵为:
\n">;
for
{
for
printf<"%-10f",R[i][j]>;
printf<"\n">;
}
for
for
Q[k][t]=H2[k][t];//此时Q为最终正交矩阵Q的转置
matrix_vec;//Q矩阵为正交阵,故Q的逆就等于Q装置,求解Q的逆
matrix_turn;//转置,求解Q
printf<"*********************************\n">;
printf<"系数矩阵A进行QR分解后的Q矩阵为:
\n">;
for
{
for
printf<"%-10f",Q[i][j]>;
printf<"\n">;
}
solution
}
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- LU QR 解法 线性方程组