北航数值分析第二次大作业QR分解.docx
《北航数值分析第二次大作业QR分解.docx》由会员分享，可在线阅读，更多相关《北航数值分析第二次大作业QR分解.docx（25页珍藏版）》请在冰豆网上搜索。
北航数值分析第二次大作业QR分解
《数值分析A》
计算实习题目二
姓名
学号
联系方式
班级
指导教师
2012年10月
数值分析第二次大作业
一、算法设计方案
首先构造矩阵A，利用Householder矩阵对矩阵A作相似变换，把A化为拟上三角矩阵A（n-1），算法如课本P61。
使用QR分解法对矩阵A（n-1）进行QR分解，算法如课本P59,进而求出所得矩阵的Q、R、RQ矩阵。
然后对A（n-1）进行带双步位移的QR分解求矩阵的全部特征值，采用以下几步进行：
第一步：
判断是否am,m-1（k）ε,若不是,则进入第四步。
若是，则am,m-1（k）为特征值，m=m-1，若m=1，则进入第二步，若m=2进入第三步，否则转第四步。
第二步：
m=1，则a11（k）为特征值，转向结束步。
第三步：
m=2,则可以求出A的两个特征值s1和s2，转向结束步。
第四步：
判断是否am-1,m-2（k）ε，若不是，进入第五步。
若是，则得到A的两个特征值s1和s2，令m=m-2，若m=1，进入第二步，若m=2进入第三步，否则进入第一步。
第五步：
判断是否达到循环上限，若达到，则结束，否则进入第六步。
第六步：
对A进行双步位移QR分解，这里的算法如课本P64，分解后转向计数步。
计数步：
对循环次数进行计数，并转向第一步。
结束步：
显示所求得的特征值。
最后对实特征值利用列主元高斯消元法求解其对应的特征向量，算法如课本p17.
二、源程序代码
#include
#include
#include
inti,j,k,l,m;//定义外部变量
doubled,h,b,c,t,s;
doubleA[10][10],AA[10][10],R[10][10],Q[10][10],RQ[10][10];
doubleX[10][10],Y[10][10],Qt[10][10],M[10][10];
doubleU[10],P[10],T[10],W[10],Re[10]={0},Im[10]={0};
doubleepsilon=1e-12;
voidmain（）
{
voidQuasiuppertriangular（doubleA[][10]）;
voidQRdecomposition（doubleA[][10]）;
voidDoublestepsQR（doubleA[][10]）;
inti,j;
for（i=0;i<10;i++）
{
for（j=0;j<10;j++）
{
A[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
Q[i][j]=0;
AA[i][j]=A[i][j];
}
A[i][i]=1.5*cos（2.2*（i+1））;
AA[i][i]=A[i][i];
}
Quasiuppertriangular（A）;//调用拟上三角化函数
printf（"\nA经过拟上三角化矩阵为：
\n\n"）;
for（i=0;i<10;i++）//输出拟上三角化矩阵
{
for（j=0;j<10;j++）
{
printf（"%.12e",A[i][j]）;//输出拟上三角化矩阵
}
printf（"\n\n"）;
}
QRdecomposition（A）;//调用QR分解函数
printf（"进行QR分解后，R矩阵为：
\n\n"）;//输出R矩阵
for（i=0;i<10;i++）
{
for（j=0;j<10;j++）
{
printf（"%.12e",R[i][j]）;
}
printf（"\n\n"）;
}
printf（"Q矩阵为:
\n\n"）;//输出Q矩阵
for（i=0;i<10;i++）
{for（j=0;j<10;j++）
{printf（"%.12e",Q[i][j]）;
}
printf（"\n\n"）;
}
printf（"RQ矩阵为：
\n\n"）;//输出RQ矩阵
for（i=0;i<10;i++）
{
for（j=0;j<10;j++）
{printf（"%.12e",RQ[i][j]）;
}
printf（"\n\n"）;
}
DoublestepsQR（A）;//调用双步位移函数
printf（"\n\n特征值实部依次为：
\n\n"）;//输出特征值实部
for（j=0;j<10;j++）
{
printf（"%.12e",Re[j]）;
}
printf（"\n\n特征值虚部依次为：
\n\n"）;//输出特征值虚部
for（j=0;j<10;j++）
{
printf（"%.12e",Im[j]）;
}
//按行输出特征向量
printf（"\n\n按行输出实特征根相应特征向量为：
\n\n"）;
for（i=0;i<10;i++）
{
if（i==1||i==2||i==5||i==6）
{
continue;
}
for（j=0;j<10;j++）
{
printf（"%.12e",X[i][j]）;
}
printf（"\n\n"）;
}
getchar（）;
}
voidQuasiuppertriangular（doubleA[][10]）
{
for（j=0;j<8;j++）
{
for（i=0;i<10;i++）
{
U[i]=0;
P[i]=0;
T[i]=0;
W[i]=0;
}
m=0;
for（i=j+2;i<10;i++）
{
if（A[i][j]!
=0）
{
m=m+1;
}
}
if（m==0）
{
continue;
}
d=0;
for（i=j+1;i<10;i++）
{
d=d+pow（A[i][j],2）;
}
d=sqrt（d）;
c=-d;
if（A[j+1][j]<=0）
{
c=d;
}
h=c*（c-A[j+1][j]）;
U[j+1]=A[j+1][j]-c;
for（i=j+2;i<10;i++）
{
U[i]=A[i][j];
}
for（i=0;i<10;i++）
{
for（k=0;k<10;k++）
{
P[i]=P[i]+U[k]*A[k][i];
}
P[i]=P[i]/h;
}
t=0;
for（i=0;i<10;i++）
{
for（k=0;k<10;k++）
{
T[i]=T[i]+U[k]*A[i][k];
}
T[i]=T[i]/h;
t=t+P[i]*U[i];
}
t=t/h;
for（i=0;i<10;i++）
{
W[i]=T[i]-t*U[i];
for（k=0;k<10;k++）
{
A[i][k]=A[i][k]-W[i]*U[k]-U[i]*P[k];
if（abs（A[i][k]）<1e-12）
{
A[i][k]=0;
}
}
}
}
}
voidQRdecomposition（doubleA[][10]）
{
for（i=0;i<10;i++）
{
for（j=0;j<10;j++）
{
RQ[i][j]=0;
Q[i][j]=0;
R[i][j]=A[i][j];
}
Q[i][i]=1;
}
for（j=0;j<9;j++）
{
for（i=0;i<10;i++）
{
U[i]=0;
P[i]=0;
W[i]=0;
}
m=0;
for（i=j+1;i<10;i++）
{
if（R[i][j]!
=0）
{
m=m+1;
}
}
if（m==0）
{
continue;
}
d=0;
for（i=j;i<10;i++）
{
d=d+pow（R[i][j],2）;
}
d=sqrt（d）;
c=-d;
if（R[j][j]<=0）
{
c=d;
}
h=c*（c-R[j][j]）;
U[j]=R[j][j]-c;
for（i=j+1;i<10;i++）
{
U[i]=R[i][j];
}
for（i=0;i<10;i++）
{
for（k=0;k<10;k++）
{
W[i]=W[i]+U[k]*Q[i][k];
}
}
for（i=0;i<10;i++）
{
for（k=0;k<10;k++）
{
Q[i][k]=Q[i][k]-（（W[i]*U[k]）/h）;
}
}
for（i=0;i<10;i++）
{
for（k=0;k<10;k++）
{
P[i]=P[i]+U[k]*R[k][i];
}
P[i]=P[i]/h;
}
for（i=0;i<10;i++）
{
for（k=0;k<10;k++）
{
R[i][k]=R[i][k]-U[i]*P[k];
if（abs（R[i][k]）{
R[i][k]=0;
}
}
}
}
for（i=0;i<10;i++）//计算A（n+1）=RQ
{
for（j=0;j<10;j++）
{
for（k=0;k<10;k++）
{
RQ[i][j]=RQ[i][j]+R[i][k]*Q[k][j];
}
}
}
}
//双步位移法计算特征值特征向量函数
voidDoublestepsQR（doubleA[][10]）
{
intL=1000,m=9;//定义最大循环次数
for（i=0;i{
for（;m>-1;）
{
if（abs（A[m][m-1]）<=epsilon）
{
Re[m]=A[m][m];
m=m-1;//降阶
if（m==0）//4
{
Re[0]=A[0][0];
break;
}
if（m==-1）
{
break;
}
if（m>1）
{
continue;
}
}
b=-A[m][m]-A[m-1][m-1];//5
c=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];
if（m==1）//6
{
if（（b*b-4*c）>=0）
{
Re[m]=（-b+sqrt（b*b-4*c））/2;
Re[m-1]=（-b-sqrt（b*b-4*c））/2;
}
if（（b*b-4*c）<0）
{
Re[m]=-b/2;Im[m]=sqrt（4*c-b*b）/2;
Re[m-1]=-b/2;Im[m-1]=-sqrt（4*c-b*b）/2;
}
m=m-1;//循环出口条件
break;
}
if（（m>1）&&（abs（A[m-1][m-2]）>epsilon））//8
{
if（i==L-1）
{
printf（"Noresults!
\n"）;
m=0;//循环出口条件
break;
}
break;
}
if（（m>1）&&（abs（A[m-1][m-2]）<=epsilon））//7
{
if（（b*b-4*c）>0）
{
Re[m]=（-b+sqrt（b*b-4*c））/2;
Re[m-1]=（-b-sqrt（b*b-4*c））/2;
}
if（（b*b-4*c）<0）
{
Re[m]=-b/2;Im[m]=sqrt（4*c-b*b）/2;
Re[m-1]=-b/2;Im[m-1]=-sqrt（4*c-b*b）/2;
}
m=m-2;//降阶
if（m>0）
{
continue;
}
if（m==0）
{
Re[0]=A[0][0];
break;
}
}
}
if（m<=0）
{
break;
}
s=A[m-1][m-1]+A[m][m];//9
t=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];
for（j=0;j<10;j++）
{
for（k=0;k<10;k++）
{
Qt[j][k]=0;
Q[j][k]=0;
M[j][k]=0;
X[j][k]=0;
Y[j][k]=0;
}
}
for（j=0;j{
for（k=0;k{
for（l=0;l{
M[j][k]=M[j][k]+A[j][l]*A[l][k];
}
}
}
for（j=0;j{
for（k=0;k{
M[j][k]=M[j][k]-s*A[j][k];
}
M[j][j]=M[j][j]+t;
}
//调用QR分解函数对M矩阵进行分解并传递参数矩阵Q
QRdecomposition（M）;
for（j=0;j<10;j++）
{
for（k=0;k<10;k++）
{
Qt[j][k]=Q[k][j];
}
}
for（j=0;j{
for（k=0;k{
for（l=0;l{
X[j][k]=X[j][k]+Qt[j][l]*A[l][k];
}
}
}
for（j=0;j{
for（k=0;k{
for（l=0;l{
Y[j][k]=Y[j][k]+X[j][l]*Q[l][k];
}
}
}
for（j=0;j<10;j++）
{
for（k=0;k<10;k++）
{
A[j][k]=Y[j][k];
}
}
}
//应用列主元高斯消元法计算实部特征向量
for（l=0;l<10;l++）
{
if（l==1||l==2||l==5||l==6）
{
continue;
}
for（k=0;k<10;k++）
{
for（m=0;m<10;m++）
{
A[k][m]=AA[k][m];
}
A[k][k]=A[k][k]-Re[l];
}
for（j=0;j<9;j++）
{
m=j;
for（i=j+1;i<10;i++）
{
if（abs（A[i][j]）>abs（A[m][j]））
{
m=i;
}
}
for（k=j;k<10;k++）
{
Y[j][k]=A[j][k];
A[j][k]=A[m][k];
A[m][k]=Y[j][k];
}
for（k=j+1;k<10;k++）
{
b=A[k][j]/A[j][j];
for（i=j;i<10;i++）
{
A[k][i]=A[k][i]-A[j][i]*b;
}
}
}
X[l][9]=1;
for（i=8;i>=0;i--）
{
c=0;
for（j=i+1;j<10;j++）
{
c=c+A[i][j]*X[l][j];
}
X[l][i]=-c/A[i][i];
}
}
}
三、程序输出结果、
升级会员 