北航数值分析第二次大作业QR分解Word格式.docx
《北航数值分析第二次大作业QR分解Word格式.docx》由会员分享，可在线阅读，更多相关《北航数值分析第二次大作业QR分解Word格式.docx（25页珍藏版）》请在冰豆网上搜索。
对循环次数进行计数，并转向第一步。
结束步：
显示所求得的特征值。
最后对实特征值利用列主元高斯消元法求解其对应的特征向量，算法如课本p17.
二、源程序代码
#include<
stdio.h>
math.h>
string.h>
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<
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"
）;
i++）//输出拟上三角化矩阵
{
for（j=0;
{
printf（"
%.12e"
A[i][j]）;
//输出拟上三角化矩阵
}
printf（"
}
QRdecomposition（A）;
//调用QR分解函数
进行QR分解后，R矩阵为：
//输出R矩阵
i++）
{
for（j=0;
{
printf（"
R[i][j]）;
printf（"
Q矩阵为:
//输出Q矩阵
for（i=0;
i++）
{for（j=0;
{printf（"
Q[i][j]）;
RQ矩阵为：
//输出RQ矩阵
{printf（"
RQ[i][j]）;
}
DoublestepsQR（A）;
//调用双步位移函数
\n\n特征值实部依次为：
//输出特征值实部
j++）
{
Re[j]）;
printf（"
\n\n特征值虚部依次为：
\n\n"
//输出特征值虚部
Im[j]）;
//按行输出特征向量
\n\n按行输出实特征根相应特征向量为：
for（i=0;
if（i==1||i==2||i==5||i==6）
{
continue;
}
for（j=0;
{
printf（"
X[i][j]）;
printf（"
}
getchar（）;
}
voidQuasiuppertriangular（doubleA[][10]）
8;
U[i]=0;
P[i]=0;
T[i]=0;
W[i]=0;
m=0;
for（i=j+2;
{
if（A[i][j]!
=0）
m=m+1;
}
if（m==0）
continue;
d=0;
for（i=j+1;
d=d+pow（A[i][j],2）;
d=sqrt（d）;
c=-d;
if（A[j+1][j]<
c=d;
h=c*（c-A[j+1][j]）;
U[j+1]=A[j+1][j]-c;
U[i]=A[i][j];
for（k=0;
k<
k++）
P[i]=P[i]+U[k]*A[k][i];
P[i]=P[i]/h;
t=0;
T[i]=T[i]+U[k]*A[i][k];
T[i]=T[i]/h;
t=t+P[i]*U[i];
t=t/h;
W[i]=T[i]-t*U[i];
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]）
{
RQ[i][j]=0;
Q[i][j]=0;
R[i][j]=A[i][j];
Q[i][i]=1;
9;
if（R[i][j]!
for（i=j;
d=d+pow（R[i][j],2）;
if（R[j][j]<
h=c*（c-R[j][j]）;
U[j]=R[j][j]-c;
U[i]=R[i][j];
W[i]=W[i]+U[k]*Q[i][k];
Q[i][k]=Q[i][k]-（（W[i]*U[k]）/h）;
P[i]=P[i]+U[k]*R[k][i];
R[i][k]=R[i][k]-U[i]*P[k];
if（abs（R[i][k]）<
epsilon）
R[i][k]=0;
}
i++）//计算A（n+1）=RQ
for（k=0;
RQ[i][j]=RQ[i][j]+R[i][k]*Q[k][j];
}
//双步位移法计算特征值特征向量函数
voidDoublestepsQR（doubleA[][10]）
intL=1000,m=9;
//定义最大循环次数
L;
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）
if（m>
1）
}
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）>
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;
}
//循环出口条件
break;
if（（m>
1）&
&
（abs（A[m-1][m-2]）>
epsilon））//8
if（i==L-1）
Noresults!
\n"
m=0;
（abs（A[m-1][m-2]）<
=epsilon））//7
m=m-2;
//降阶
if（m>
continue;
if（m==0）
Re[0]=A[0][0];
if（m<
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（k=0;
Qt[j][k]=0;
Q[j][k]=0;
M[j][k]=0;
X[j][k]=0;
Y[j][k]=0;
m+1;
for（l=0;
l<
l++）
M[j][k]=M[j][k]+A[j][l]*A[l][k];
for（j=0;
M[j][k]=M[j][k]-s*A[j][k];
M[j][j]=M[j][j]+t;
//调用QR分解函数对M矩阵进行分解并传递参数矩阵Q
QRdecomposition（M）;
Qt[j][k]=Q[k][j];
X[j][k]=X[j][k]+Qt[j][l]*A[l][k];
Y[j][k]=Y[j][k]+X[j][l]*Q[l][k];
A[j][k]=Y[j][k];
//应用列主元高斯消元法计算实部特征向量
for（l=0;
l++）
if（l==1||l==2||l==5||l==6）
continue;
for（m=0;
m<
m++）
A[k][m]=AA[k][m];
A[k][k]=A[k][k]-Re[l];
m=j;
for（i=j+1;
if（abs（A[i][j]）>
abs（A[m][j]））
m=i;
}
for（k=j;
Y[j][k]=A[j][k];
A[j][k]=A[m][k];
A[m][k]=Y[j][k];
for（k=j+1;
k++）
b=A[k][j]/A[j][j];
for（i=j;
{
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;
c=c+A[i][j]*X[l][j];
X[l][i]=-c/A[i][i];
三、程序输出结果、
升级会员 