线性代数方程组求解.docx
《线性代数方程组求解.docx》由会员分享，可在线阅读，更多相关《线性代数方程组求解.docx（21页珍藏版）》请在冰点文库上搜索。
线性代数方程组求解
线性代数方程组求解
一、实验要求
编程求解方程组:
方程组１:
方程组2:
方程组３:
要求:
用Ｃ/C++语言实现如下函数:
1.ｂooｌlu（dｏublｅ*ａ, int*ｐiｖot,intn）;
实现矩阵的LU分解。
ｐiｖｏｔ为输出参数,ｐｉvoｔ［0,n）中存放主元的位置排列。
函数成功时返回false,否则返回trｕe。
2.bｏoｌguasｓ（doｕblｅcoｎst*ｌu,inｔconst*p,dｏuble＊ b,intn）;
求线代数方程组的解
设矩阵Lunxｎ为某个矩阵anｘｎ的LＵ分解,在内存中按行优先次序存放。
p[0,n）为ＬU分解的主元排列。
b为方程组Ax=ｂ的右端向量。
此函数计算方程组Ax=b的解,并将结果存放在数组b[0,n）中。
函数成功时返回falｓｅ,否则返回trｕe。
3、voｉｄ qｒ（dｏuble*a,ｄoubｌe*d,intn）;矩阵的QR分解
假设数组aｎxn在内存中按行优先次序存放。
此函数使用HｏuｓeHｏlｄer变换将其就地进行QR分解。
d为输出参数,d[0,ｎ）中存放QR分解的上三角对角线元素。
４、booｌhｓｈlｄ（dｏubleconｓt*ｑr, doublｅcoｎsｔ*d,ｄoubｌe*b, intｎ）;求线代数方程组的解
设矩阵ｑrnｘn为某个矩阵anxn的QＲ分解,在内存中按行优先次序存放。
d [0,n）为QR分解的上三角对角线元素。
b为方程组Ax＝b的右端向量。
函数计算方程组Ａx=b的解,并将结果存放在数组b[0,n）中。
函数成功时返回falｓe,否则返回tｒue。
二、问题分析
求解线性方程组Ａx=ｂ,其实质就就是把它的系数矩阵A通过各种变换成一个下三角或上三角矩阵,从而简化方程组的求解。
因此,在求解线性方程组的过程中,把系数矩阵A变换成上三角或下三角矩阵显得尤为重要,然而矩阵A的变换通常有两种分解方法:
LU分解法与QR分解法。
1、LU分解法:
将A分解为一个下三角矩阵L与一个上三角矩阵U,即:
A＝LＵ,
其中 L=
U=
2、ＱＲ分解法:
将A分解为一个正交矩阵Q与一个上三角矩阵Ｒ,即:
Ａ=ＱR
三、实验原理
解Ax=b的问题就等价于要求解两个三角形方程组:
ﻫ ⑴Ly＝b,求y;
⑵Ｕx=y,求x、
设A为非奇异矩阵,且有分解式A＝LU, Ｌ为单位下三角阵,U为上三角阵。
Ｌ,U的元素可以有n步直接计算定出。
用直接三角分解法解Ａx＝b（要求A的所有顺序主子式都不为零）的计算公式:
①
i=２,３,…,n.
计算U的第r行,L的第r列元素（i=2,3,…,n）:
ﻫ②
, i=r,ｒ+1,…,n;
③
i=r+1,…,n,且r≠n、ﻫ求解Ｌy=b,Ux=y的计算公式;
④ﻩ
⑤
四、实验步骤
1>将矩阵A保存进计算机中,再定义２个空矩阵L,Ｕ以便保存求出的三角矩阵的值。
利用公式①,②,③将矩阵A分解为LＵ,Ｌ为单位下三角阵,Ｕ为上三角阵。
2>可知计算方法有三层循环。
先通过公式①计算出U矩阵的第一行元素
与L矩阵的第一列元素
。
再根据公式②与③,与上次的出的值,求出矩阵其余的元素,每次都要三次循环,求下一个元素需要上一个结果。
3>先由公式④,Ly=ｂ
求出y,因为L为下三角矩阵,所以由第一行开始求y、
4>再由公式⑤,Ux=y
求出x,因为Ｕ为上三角矩阵,所以由最后一行开始求x、
五、程序流程图
1、LＵ分解法
2、ＱR分解法
ﻬ六、实验结果
1、LU分解法
方程组１:
方程组2:
方程组3:
2、QR分解法
方程组1:
方程组2:
ﻩ方程组3:
七、实验总结
为了求解线性方程组,我们通常需要一定的解法。
其中一种解法就就是通过矩阵的三角分解来实现的,属于求解线性方程组的直接法。
在不考虑舍入误差下,直接法可以用有限的运算得到精确解,因此主要适用于求解中小型稠密的线性方程组。
1、三角分解法
三角分解法就是将Ａ矩阵分解成一个上三角形矩阵U与一个下三角形矩阵L,这样的分解法又称为LU分解法。
它的用途主要在简化一个大矩阵的行列式值的计算过程,求反矩阵与求解联立方程组。
不过要注意这种分解法所得到的上下三角形矩阵并非唯一,还可找到数个不同 的一对上下三角形矩阵,此两三角形矩阵相乘也会得到原矩阵。
2、 QR分解法
QＲ分解法就是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为QＲ分解法。
在编写这两个程序过程中,起初遇到不少麻烦!
虽然课上老师反复重复着:
“算法不难的,It＇sｓoｅasy！
”但就是当自己实际操作时,感觉并不就是那么容易。
毕竟就是要把实际的数学问题转化为计算机能够识别的编程算法,所以在编写程序之前我们仔细认真的把所求解的问题逐一进行详细的分析,最终转化为程序段。
每当遇到问题时,大家或许有些郁闷,但最终还就是静下心来反复仔细的琢磨,一一排除了错误,最终完成了本次实验。
回头一想原来编个程序其实也没有想象的那么复杂,只要思路清晰,逐步分析,就可以慢慢搞定了。
附源代码:
＃ｉnclude
#inclｕｄe＜stｄｉo、h>
#inｃlude<ｍath.h>
＃incｌudｅ<ｓtdlｉb、h>
#iｎclｕde
usｉng nａmespaｃestd;
bｏollu（dｏuble*a,ｉｎt＊ｐivot,intｎ）;/／矩阵LU分解
boｏｌｇuａss（ｄoubleconｓt*lu,intｃｏnst＊ p, doubｌe*b,intn）;/／求线性代数方程组的解
ｖoｉd qｒ（double*a,double* ｄ, intn）; //矩阵的QR分解
ｂoｏｌ hｏuseholｄer（doublｅcoｎst*qｒ,doｕｂlｅconsｔ*d, doubｌｅ*ｂ, inｔ ｎ）;
inｔｍain（）
{
iｎt n=0;
ｉnｔtｅmp=0;
ｂooｌfｌａg＝falsｅ;
ｄoｕbleexpct=0;/／误差期望值
dｏuｂｌe devsｑ＝0;//误差的方差
int*P=NULL;
dｏuble* D＝NＵLL;
doｕｂleＡ[］={1,1/2、0,１/3、0, １/4、0,１/５、0,1/６.０,
１／2、０,１/３、0,１／4、0, 1/5、0, １／６、０,1/7、0,
1/3、０,1／4.0,1/5.0,１/6、0,1/７.0,１/８、0,
1/4、0,1/5、０,1/6.０,1/７、0,1/8、0,1／9、0,
1／５、０,1／6、0,1/7、0,1/８、0, 1／９、0,1/１0、０,
１/6、0,1/７.0,1/8、0, 1/9、０,１／10、0,1／11、０
};
double B［]=｛ 1＋ 1/2、0+1／3、０+1/4、０+1/5、0+1/6、0,
1/2.0＋1/3、0+1／４.0＋1／５.０+1/6、0+１/７、０,
１/3、0+1/4、0+1／5、0+ 1/６.０+1/7.０+1/8、０,
1／4.0+１/5、0+1/6、0+1／７、0+1/8.0+１/9、0,
1/5、0+１/6、0+ 1／７、0+1/8、０+1／9、0＋1/1０.0,
１/6、０＋1/7、０+1／８、0+1／9、0＋1/10.０+1／11、0
｝;
n=6;
P= （iｎｔ*）ｍallｏc（sｉzeｏf（int）*n）;
D=（ｄouｂle*）mａｌloc（sｉzeof（double）*n）;
for （intｉ=０;i<ｎ;i+＋）
{
P[i］=D[ｉ]=０;
}
cｏｕt<<"Please ｃｈoｏseｔhｅｗaytoｓolｖe（0:
lu,1:
qr）:
";
ciｎ>>flaｇ;
if（!
flag）
{
coｕｔ＜<"矩阵ＬU分解:
\n＂;
lu（A,P,n）;//矩阵LU分解
ｆor（ｉnti=0; i {
for（intj=0; ｊ prinｔf（"%ｆ\ｔ",A［i]）;
cout< }
cｏut<<＂矩阵LＵ分解的主元位置:
\n";
for（inti=0;i cout＜<Ｐ［i]＜<"＼t";
ｃout<＜"＼ｎGｕａss求线性代数方程组求解:
\n";
gｕaｓs（A,P,B,n）;/／求线性代数方程组的解
ｆｏr（inti=0; i＜n;i++）
{
ｉf（ｉ=＝3）
ｃout＜＜ｅndl;
ｐrintf（"%、16f＼ｔ＂,Ｂ[i］）;
}
coｕt<<"\ｎGuass解法的误差为:
\n";
for（iｎtｉ=0;ｉ＜n; i+＋）
｛
ｉf（ｉ==３）
couｔ＜<ｅndｌ;
pｒｉnｔｆ（＂%、1６ｆ\t",B［i]－1）;
expct=expcｔ+Ｂ［i］-1;
｝
prｉntf（"＼n误差期望值为%、16f\t\n",expct／6）;
cout＜<ｅndl;
}
elｓe
{
ｃｏｕt<<"矩阵ＱＲ分解:
＼n";
qr（Ａ,Ｄ,n）;//矩阵qr分解
fｏr（inti=０; i ｛
ｆoｒ（ｉｎｔ j=0;j<ｎ;j＋+）
prinｔｆ（＂%ｄ\t",A［ｉ]）;
cout<<ｅｎdl;
}
ｃout<<"QR分解的上三角矩阵对角线元素:
＼ｎ";
for（inti=0;i pｒintf（＂%f\t",D[i］）;
ｃoｕt<<"\nＨｏｕｓｅholder线性代数方程组求解:
\n";
houseｈolder（A,Ｄ,Ｂ,ｎ）;//求线性代数方程组的解
fｏｒ（inｔi=0;i {
if（i=＝3）
cｏut＜<ｅndl;
printｆ（"%.16ｆ\ｔ",B[i]）;
}
cout＜<"\nHousｅholdｅr解法的误差为:
\n＂;
foｒ（ｉntｉ=０;i<ｎ;i＋+）
{
if（i＝=３）
coｕt＜ ｐrintf（"％、16ｆ\t＂,B[i］－1）;
expcｔ=expcｔ+B[ｉ]-１;
}
printf（"\ｎ误差期望值为%、１6f\ｔ\n",exｐｃt／6）;
cout< }
ｒeturn0;
}
boｏｌlu（doｕbｌe*a,iｎt* pivot, int ｎ）//矩阵LU分解
｛
iｎti,j,ｋ;
doｕbｌｅｍａx,temp;
ｍａｘ =0;
ｔｅmp=0;
for（ｉ=0;i {
//选出i列的主元,记录主元位置
maｘ=fabs（ａ[ｎ*ｉ＋ ｉ]）;
pｉvｏｔ[ｉ]=i;
for（j=i+１; j {
if（fabs（a[n＊j＋i]）>max）
{
ｍaｘ=fabs（a［n*j ＋ ｉ］）;
piｖot［ｉ]=j;
}
}
//对第i列进行行变换,使得主元在对角线上
if（pivｏt[i]!
＝i）
{
fｏr（j=i;ｊ {
teｍｐ=ａ[n＊i＋ｊ];
a［n＊i+j]=a[n*ｐｉvoｔ[ｉ]＋j];
a[n*ｐｉvoｔ[i]+ｊ]=ｔｅmp;
}
}
ｆor（j=ｉ+1;ｊ a[n＊j +ｉ]=a[n*j+i]/ａ[n＊i+ｉ];
for（ｊ=i+1;j<ｎ;j＋+）/／计算上三角U
foｒ（k=i+1;k a[n＊j+k]=a［n*j＋k]－a[n*j+i］*a[n＊i+k];
}
／/计算下三角L
for（i=0; i<ｎ-2;i＋+）/／ｉ行k列
foｒ（ｋ=i＋1; k ｛
temｐ＝ａ[n*pivot[k］+i］;
a[n*ｐivｏt[ｋ]+i］=a［k*n+i];
a［k*n+ i]=teｍｐ;
}
rｅtuｒnfａlｓe;
}
booｌ guasｓ（ｄoubleconst＊lu,intconst*ｐ,dｏuｂｌｅ*b,ｉnt n）//求线性代数方程组的解
{
intｉ,j,k;
douｂle ｔemp;
/／按ｑivot对ｂ行变换,与LU匹配
ｆor（i=０;ｉ ｛
teｍp=b[ｐ[i]];
ｂ[p[i]］=b[i];
b[i]=temp;
}
//Ｌy=b,将y的内容放入b
ｆor（ｉ=0;i ｆor（j＝0; j b[i]＝b[i]－lｕ[n*i+j]＊b[j];
/／Uy=x,将ｘ的内容放入ｂ
for（ｉ＝n－1;ｉ＞=0;i-－）
｛
fｏｒ（j=n－１; ｊ>i;ｊ--）
ｂ[i]＝b[i]-lu[n＊ｉ+ｊ］＊b[j];
ｂ[i］=b［ｉ]／lu[n*ｉ+i];
}
retuｒｎfaｌsｅ;
}
voidqr（doubｌｅ*a,dｏuｂle*d,int n）//矩阵的ＱR分解
{
iｎｔｉ,ｊ,l,ｋ;
doubletem,ｍ;
douｂlｅ＊ｔemp;
tｅmｐ＝（double*）ｍalloc（sｉzｅｏf（doublｅ）*n）;
ｆoｒ（i=0; i{
/／获得tｅm值
m＝0;
for（j＝i;j ｍ=m+a[ｎ*j+ i]*a[n*ｊ+ｉ］;
ｉf（a［ｎ*i＋i ]>0）
m=-sqrt（ｍ）;
elｓe
m=ｓqrt（ｍ）;
//获得temｐ放入矩阵,并存主元ｄ
tem =０ ;
d［i]＝m;
ａ[n*i＋i] = ａ[n＊ｉ+ｉ］－ｍ;
foｒ（j＝i;ｊ<=n-１;j+＋）
tem=tem+ a[ｎ*j+i］*a［n*ｊ+i];
teｍ=sqｒｔ（ｔem）;
ｆor（ｊ＝i;j＜=ｎ-1;j+＋）
a[ｎ*j +i]= a[n*j＋i］/tｅm;
//调整矩阵
ｆoｒ（k＝i+1;k {
fｏｒ（j=i;j {
ｔｅm=0 ;
for（l＝ｉ; ｌ tem=tｅm+ a[ｎ*j+ i]*a[n＊l+i]*a[n*l＋ｋ];
tｅｍp[j]= a[j＊ｎ＋k] － ２＊tem;
}
foｒ（j=i;j a［j*n＋k]＝ temｐ［j];
｝
｝
d[n－１］ =a[（n-１）*n+n－1］;
}
ｂｏｏlhousｅholder（ｄoubｌecoｎst*qr, double cｏnst*d, ｄouｂle＊ｂ, int n）//求线性代数方程组的解
｛
int ｉ,j,l;
ｄouｂleｒem;
douｂｌe＊tｅmp;
temp＝（ｄoｕbｌe＊）maｌｌoc（sｉｚeｏf（ｄoｕble）*n）;
for（i=０; ｉ ｛
fｏr（j＝i;j {
reｍ=０;
for（l=i;ｌ rｅｍ=rｅm ＋ qr[l*n＋i]*qｒ[j*ｎ+ｉ]*b［l];
ｔｅｍｐ［ｊ]=b[j]-2＊rem;
}
foｒ（j=ｉ; j＜n;j＋＋）
b[j］=tｅmp[j];
｝
ｆor（j=ｎ-1; j>-１; j－－）
{
fｏr（l=n-1;l！
=j;--l）
b[ｊ］=b[j] -b[l]*qr[ｊ*n＋l］;
ｂ[j]=b[j］／d［j];
｝
ｒｅtuｒn false;
}
升级会员 