线性代数方程组求解.docx
《线性代数方程组求解.docx》由会员分享，可在线阅读，更多相关《线性代数方程组求解.docx（19页珍藏版）》请在冰豆网上搜索。
线性代数方程组求解
线性代数方程组求解
一、实验要求
编程求解方程组：
方程组1：
方程组2：
方程组3：
要求：
用C/C++语言实现如下函数：
1.boollu（double*a,int*pivot,intn）;
实现矩阵的LU分解。
pivot为输出参数，pivot[0,n）中存放主元的位置排列。
函数成功时返回false，否则返回true。
2.boolguass（doubleconst*lu,intconst*p,double*b,intn）;
求线代数方程组的解
设矩阵Lunxn为某个矩阵anxn的LU分解，在内存中按行优先次序存放。
p[0,n）为LU分解的主元排列。
b为方程组Ax=b的右端向量。
此函数计算方程组Ax=b的解，并将结果存放在数组b[0,n）中。
函数成功时返回false，否则返回true。
3.voidqr（double*a,double*d,intn）;矩阵的QR分解
假设数组anxn在内存中按行优先次序存放。
此函数使用HouseHolder变换将其就地进行QR分解。
d为输出参数，d[0,n）中存放QR分解的上三角对角线元素。
4.boolhshld（doubleconst*qr,doubleconst*d,double*b,intn）;求线代数方程组的解
设矩阵qrnxn为某个矩阵anxn的QR分解，在内存中按行优先次序存放。
d[0,n）为QR分解的上三角对角线元素。
b为方程组Ax=b的右端向量。
函数计算方程组Ax=b的解，并将结果存放在数组b[0,n）中。
函数成功时返回false，否则返回true。
二、问题分析
求解线性方程组Ax=b，其实质就是把它的系数矩阵A通过各种变换成一个下三角或上三角矩阵，从而简化方程组的求解。
因此，在求解线性方程组的过程中，把系数矩阵A变换成上三角或下三角矩阵显得尤为重要，然而矩阵A的变换通常有两种分解方法：
LU分解法和QR分解法。
1、LU分解法：
将A分解为一个下三角矩阵L和一个上三角矩阵U，即：
A=LU,
其中L=
，U=
2、QR分解法：
将A分解为一个正交矩阵Q和一个上三角矩阵R,即：
A=QR
三、实验原理
解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=2,3,…,n）：
②
，i=r,r+1,…,n;
③
i=r+1,…,n,且r≠n.
求解Ly=b，Ux=y的计算公式；
④
⑤
四、实验步骤
1>将矩阵A保存进计算机中，再定义2个空矩阵L，U以便保存求出的三角矩阵的值。
利用公式①，②，③将矩阵A分解为LU，L为单位下三角阵,U为上三角阵。
2>可知计算方法有三层循环。
先通过公式①计算出U矩阵的第一行元素
和L矩阵的第一列元素
。
再根据公式②和③，和上次的出的值，求出矩阵其余的元素，每次都要三次循环，求下一个元素需要上一个结果。
3>先由公式④，Ly=b
求出y，因为L为下三角矩阵，所以由第一行开始求y.
4>再由公式⑤，Ux=y
求出x,因为U为上三角矩阵，所以由最后一行开始求x.
五、程序流程图
1、LU分解法
2、QR分解法
六、实验结果
1、LU分解法
方程组1：
方程组2：
方程组3：
2、QR分解法
方程组1：
方程组2：
方程组3：
七、实验总结
为了求解线性方程组，我们通常需要一定的解法。
其中一种解法就是通过矩阵的三角分解来实现的，属于求解线性方程组的直接法。
在不考虑舍入误差下，直接法可以用有限的运算得到精确解，因此主要适用于求解中小型稠密的线性方程组。
1、三角分解法
三角分解法是将A矩阵分解成一个上三角形矩阵U和一个下三角形矩阵L，这样的分解法又称为LU分解法。
它的用途主要在简化一个大矩阵的行列式值的计算过程，求反矩阵和求解联立方程组。
不过要注意这种分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同的一对上下三角形矩阵，此两三角形矩阵相乘也会得到原矩阵。
2、QR分解法
QR分解法是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为QR分解法。
在编写这两个程序过程中，起初遇到不少麻烦！
虽然课上老师反复重复着：
“算法不难的，It'ssoeasy!
”但是当自己实际操作时，感觉并不是那么容易。
毕竟是要把实际的数学问题转化为计算机能够识别的编程算法，所以在编写程序之前我们仔细认真的把所求解的问题逐一进行详细的分析，最终转化为程序段。
每当遇到问题时，大家或许有些郁闷，但最终还是静下心来反复仔细的琢磨，一一排除了错误，最终完成了本次实验。
回头一想原来编个程序其实也没有想象的那么复杂，只要思路清晰，逐步分析，就可以慢慢搞定了。
附源代码:
#include
#include
#include
#include
#include
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;
doubleexpct=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{
P[i]=D[i]=0;
}
cout<<"Pleasechoosethewaytosolve（0:
lu,1:
qr）:
";
cin>>flag;
if（!
flag）
{
cout<<"矩阵LU分解:
\n";
lu（A,P,n）;//矩阵LU分解
for（inti=0;i{
for（intj=0;jprintf（"%f\t",A[i]）;
cout<}
cout<<"矩阵LU分解的主元位置:
\n";
for（inti=0;icout<cout<<"\nGuass求线性代数方程组求解:
\n";
guass（A,P,B,n）;//求线性代数方程组的解
for（inti=0;i{
if（i==3）
cout<printf（"%.16f\t",B[i]）;
}
cout<<"\nGuass解法的误差为:
\n";
for（inti=0;i{
if（i==3）
cout<printf（"%.16f\t",B[i]-1）;
expct=expct+B[i]-1;
}
printf（"\n误差期望值为%.16f\t\n",expct/6）;
cout<}
else
{
cout<<"矩阵QR分解:
\n";
qr（A,D,n）;//矩阵qr分解
for（inti=0;i{
for（intj=0;jprintf（"%d\t",A[i]）;
cout<}
cout<<"QR分解的上三角矩阵对角线元素:
\n";
for（inti=0;iprintf（"%f\t",D[i]）;
cout<<"\nHouseholder线性代数方程组求解:
\n";
householder（A,D,B,n）;//求线性代数方程组的解
for（inti=0;i{
if（i==3）
cout<printf（"%.16f\t",B[i]）;
}
cout<<"\nHouseholder解法的误差为:
\n";
for（inti=0;i{
if（i==3）
cout<printf（"%.16f\t",B[i]-1）;
expct=expct+B[i]-1;
}
printf（"\n误差期望值为%.16f\t\n",expct/6）;
cout<}
return0;
}
boollu（double*a,int*pivot,intn）//矩阵LU分解
{
inti,j,k;
doublemax,temp;
max=0;
temp=0;
for（i=0;i{
//选出i列的主元,记录主元位置
max=fabs（a[n*i+i]）;
pivot[i]=i;
for（j=i+1;j{
if（fabs（a[n*j+i]）>max）
{
max=fabs（a[n*j+i]）;
pivot[i]=j;
}
}
//对第i列进行行变换，使得主元在对角线上
if（pivot[i]!
=i）
{
for（j=i;j{
temp=a[n*i+j];
a[n*i+j]=a[n*pivot[i]+j];
a[n*pivot[i]+j]=temp;
}
}
for（j=i+1;ja[n*j+i]=a[n*j+i]/a[n*i+i];
for（j=i+1;jfor（k=i+1;ka[n*j+k]=a[n*j+k]-a[n*j+i]*a[n*i+k];
}
//计算下三角L
for（i=0;ifor（k=i+1;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）//求线性代数方程组的解
{
inti,j,k;
doubletemp;
//按qivot对b行变换，与LU匹配
for（i=0;i{
temp=b[p[i]];
b[p[i]]=b[i];
b[i]=temp;
}
//Ly=b,将y的内容放入b
for（i=0;ifor（j=0;jb[i]=b[i]-lu[n*i+j]*b[j];
//Uy=x,将x的内容放入b
for（i=n-1;i>=0;i--）
{
for（j=n-1;j>i;j--）
b[i]=b[i]-lu[n*i+j]*b[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）;
for（i=0;i{
//获得tem值
m=0;
for（j=i;jm=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;
for（j=i;j<=n-1;j++）
tem=tem+a[n*j+i]*a[n*j+i];
tem=sqrt（tem）;
for（j=i;j<=n-1;j++）
a[n*j+i]=a[n*j+i]/tem;
//调整矩阵
for（k=i+1;k{
for（j=i;j{
tem=0;
for（l=i;ltem=tem+a[n*j+i]*a[n*l+i]*a[n*l+k];
temp[j]=a[j*n+k]-2*tem;
}
for（j=i;ja[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;
double*temp;
temp=（double*）malloc（sizeof（double）*n）;
for（i=0;i{
for（j=i;j{
rem=0;
for（l=i;lrem=rem+qr[l*n+i]*qr[j*n+i]*b[l];
temp[j]=b[j]-2*rem;
}
for（j=i;jb[j]=temp[j];
}
for（j=n-1;j>-1;j--）
{
for（l=n-1;l!
=j;--l）
b[j]=b[j]-b[l]*qr[j*n+l];
b[j]=b[j]/d[j];
}
returnfalse;
}
升级会员 