线性代数方程组求解文档格式.docx
《线性代数方程组求解文档格式.docx》由会员分享，可在线阅读，更多相关《线性代数方程组求解文档格式.docx（12页珍藏版）》请在冰豆网上搜索。
LU分解法和QR分解法。
1、LU分解法：
将A分解为一个下三角矩阵L和一个上三角矩阵U，即：
A=LU,
其中L=
，U=
2、QR分解法：
将A分解为一个正交矩阵Q和一个上三角矩阵R,即：
A=QR
三、实验原理
解A*=b的问题就等价于要求解两个三角形方程组：
⑴Ly=b,求y;
⑵U*=y,求*.
设A为非奇异矩阵,且有分解式A=LU,L为单位下三角阵,U为上三角阵。
L,U的元素可以有n步直接计算定出。
用直接三角分解法解A*=b〔要求A的所有顺序主子式都不为零〕的计算公式：
①
i=2,3,…，n.
计算U的第r行，L的第r列元素（i=2,3,…,n）：
②
，i=r,r+1,…,n;
③
i=r+1,…,n,且r≠n.
求解Ly=b，U*=y的计算公式；
④
⑤
四、实验步骤
1>
将矩阵A保存进计算机中，再定义2个空矩阵L，U以便保存求出的三角矩阵的值。
利用公式①，②，③将矩阵A分解为LU，L为单位下三角阵,U为上三角阵。
2>
可知计算方法有三层循环。
先通过公式①计算出U矩阵的第一行元素
和L矩阵的第一列元素
。
再根据公式②和③，和上次的出的值，求出矩阵其余的元素，每次都要三次循环，求下一个元素需要上一个结果。
3>
先由公式④，Ly=b
求出y，因为L为下三角矩阵，所以由第一行开场求y.
4>
再由公式⑤，U*=y
求出*,因为U为上三角矩阵，所以由最后一行开场求*.
五、程序流程图
1、LU分解法
2、QR分解法
六、实验结果
方程组1：
方程组2：
方程组1：
方程组3：
七、实验总结
为了求解线性方程组，我们通常需要一定的解法。
其中一种解法就是通过矩阵的三角分解来实现的，属于求解线性方程组的直接法。
在不考虑舍入误差下，直接法可以用有限的运算得到准确解，因此主要适用于求解中小型稠密的线性方程组。
1、三角分解法
三角分解法是将A矩阵分解成一个上三角形矩阵U和一个下三角形矩阵L，这样的分解法又称为LU分解法。
它的用途主要在简化一个大矩阵的行列式值的计算过程，求反矩阵和求解联立方程组。
不过要注意这种分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同的一对上下三角形矩阵，此两三角形矩阵相乘也会得到原矩阵。
2、QR分解法
QR分解法是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为QR分解法。
在编写这两个程序过程中，起初遇到不少麻烦！
虽然课上教师反复重复着：
“算法不难的，It'
ssoeasy!
〞但是当自己实际操作时，感觉并不是则容易。
毕竟是要把实际的数学问题转化为计算机能够识别的编程算法，所以在编写程序之前我们仔细认真的把所求解的问题逐一进展详细的分析，最终转化为程序段。
每当遇到问题时，大家或许有些郁闷，但最终还是静下心来反复仔细的琢磨，一一排除了错误，最终完成了本次实验。
回头一想原来编个程序其实也没有想象的则复杂，只要思路清晰，逐步分析，就可以慢慢搞定了。
附源代码:
*include<
iostream>
stdio.h>
math.h>
stdlib.h>
conio.h>
usingnamespacestd;
boollu（double*a,int*pivot,intn）;
//矩阵LU分解
boolguass（doubleconst*lu,intconst*p,double*b,intn）;
//求线性代数方程组的解
voidqr（double*a,double*d,intn）;
//矩阵的QR分解
boolhouseholder（doubleconst*qr,doubleconst*d,double*b,intn）;
intmain（）
{
intn=0;
inttemp=0;
boolflag=false;
doublee*pct=0;
//误差期望值
doubledevsq=0;
//误差的方差
int*P=NULL;
double*D=NULL;
doubleA[]={1,1/2.0,1/3.0,1/4.0,1/5.0,1/6.0,
1/2.0,1/3.0,1/4.0,1/5.0,1/6.0,1/7.0,
1/3.0,1/4.0,1/5.0,1/6.0,1/7.0,1/8.0,
1/4.0,1/5.0,1/6.0,1/7.0,1/8.0,1/9.0,
1/5.0,1/6.0,1/7.0,1/8.0,1/9.0,1/10.0,
1/6.0,1/7.0,1/8.0,1/9.0,1/10.0,1/11.0
};
doubleB[]={1+1/2.0+1/3.0+1/4.0+1/5.0+1/6.0,
1/2.0+1/3.0+1/4.0+1/5.0+1/6.0+1/7.0,
1/3.0+1/4.0+1/5.0+1/6.0+1/7.0+1/8.0,
1/4.0+1/5.0+1/6.0+1/7.0+1/8.0+1/9.0,
1/5.0+1/6.0+1/7.0+1/8.0+1/9.0+1/10.0,
1/6.0+1/7.0+1/8.0+1/9.0+1/10.0+1/11.0
n=6;
P=（int*）malloc（sizeof（int）*n）;
D=（double*）malloc（sizeof（double）*n）;
for（inti=0;
i<
n;
i++）
{
P[i]=D[i]=0;
}
cout<
<
"
Pleasechoosethewaytosolve（0:
lu,1:
qr）:
;
cin>
>
flag;
if（!
flag）
矩阵LU分解:
\n"
lu（A,P,n）;
//矩阵LU分解
for（intj=0;
j<
j++）
printf（"
%f\t"
A[i]）;
endl;
矩阵LU分解的主元位置:
P[i]<
\t"
\nGuass求线性代数方程组求解:
\n"
guass（A,P,B,n）;
if（i==3）
printf（"
%.16f\t"
B[i]）;
}
\nGuass解法的误差为:
B[i]-1）;
e*pct=e*pct+B[i]-1;
\n误差期望值为%.16f\t\n"
e*pct/6）;
else
矩阵QR分解:
qr（A,D,n）;
//矩阵qr分解
%d\t"
QR分解的上三角矩阵对角线元素:
D[i]）;
\nHouseholder线性代数方程组求解:
householder（A,D,B,n）;
//求线性代数方程组的解
\nHouseholder解法的误差为:
return0;
}
boollu（double*a,int*pivot,intn）//矩阵LU分解
inti,j,k;
doublema*,temp;
ma*=0;
temp=0;
for（i=0;
n-1;
i++）//依次对第i列进展处理
//选出i列的主元,记录主元位置
ma*=fabs（a[n*i+i]）;
pivot[i]=i;
for（j=i+1;
j++）//j表示第j行
if（fabs（a[n*j+i]）>
ma*）
ma*=fabs（a[n*j+i]）;
pivot[i]=j;
//对第i列进展行变换，使得主元在对角线上
if（pivot[i]!
=i）
for（j=i;
j++）//ij与pivot[i]j换只用对上三角进展处理
temp=a[n*i+j];
a[n*i+j]=a[n*pivot[i]+j];
a[n*pivot[i]+j]=temp;
j++）//Pi局部下三角L
a[n*j+i]=a[n*j+i]/a[n*i+i];
j++）//计算上三角U
for（k=i+1;
k<
k++）
a[n*j+k]=a[n*j+k]-a[n*j+i]*a[n*i+k];
//计算下三角L
for（i=0;
n-2;
i++）//i行k列
k++）
temp=a[n*pivot[k]+i];
a[n*pivot[k]+i]=a[k*n+i];
a[k*n+i]=temp;
returnfalse;
}
boolguass（doubleconst*lu,intconst*p,double*b,intn）//求线性代数方程组的解
doubletemp;
//按qivot对b行变换，与LU匹配
i++）//貌似错误在这里哦
temp=b[p[i]];
b[p[i]]=b[i];
b[i]=temp;
//Ly=b,将y的内容放入b
i++）
for（j=0;
i;
b[i]=b[i]-lu[n*i+j]*b[j];
//Uy=*,将*的内容放入b
for（i=n-1;
i>
=0;
i--）
for（j=n-1;
j>
j--）
b[i]=b[i]/lu[n*i+i];
returnfalse;
voidqr（double*a,double*d,intn）//矩阵的QR分解
inti,j,l,k;
doubletem,m;
double*temp;
temp=（double*）malloc（sizeof（double）*n）;
//获得tem值
m=0;
m=m+a[n*j+i]*a[n*j+i];
if（a[n*i+i]>
0）
m=-sqrt（m）;
else
m=sqrt（m）;
//获得temp放入矩阵，并存主元d
tem=0;
d[i]=m;
a[n*i+i]=a[n*i+i]-m;
=n-1;
tem=tem+a[n*j+i]*a[n*j+i];
tem=sqrt（tem）;
a[n*j+i]=a[n*j+i]/tem;
//调整矩阵
k<
{
for（l=i;
l<
l++）
tem=tem+a[n*j+i]*a[n*l+i]*a[n*l+k];
temp[j]=a[j*n+k]-2*tem;
a[j*n+k]=temp[j];
d[n-1]=a[（n-1）*n+n-1];
boolhouseholder（doubleconst*qr,doubleconst*d,double*b,intn）//求线性代数方程组的解
inti,j,l;
doublerem;
rem=0;
l<
rem=rem+qr[l*n+i]*qr[j*n+i]*b[l];
temp[j]=b[j]-2*rem;
b[j]=temp[j];
-1;
for（l=n-1;
l!
=j;
--l）
b[j]=b[j]-b[l]*qr[j*n+l];
b[j]=b[j]/d[j];
升级会员 