北航研究生数值分析大作业二.docx
《北航研究生数值分析大作业二.docx》由会员分享，可在线阅读，更多相关《北航研究生数值分析大作业二.docx（22页珍藏版）》请在冰豆网上搜索。
北航研究生数值分析大作业二
数值分析
—计算实习作业二
学院：
17系
专业：
精密仪器及机械
姓名：
张大军
学号：
DY1417114
一、程序设计方案
程序设计方案流程图如图1所示。
（注：
由本人独立完成，并且有几处算法很巧妙）
2、程序源代码
#include
#include
#include
#defineN10
#defineE1.0e-12
#defineMAX10000
intmain（）
{
voidnishangsanjiaohua（double（*A）[10]）;
voidQRfenjie（double（*A）[10],double（*Q）[N],double（*R）[N]）;
voidzhengli（double（*A）[10]）;
voidsubuQR（double（*A）[10],double*RR,double*II）;
voidtezhengxl（double（*a）[N],doubleT）;
doubleA[10][10]={0},Q[10][10]={0},R[10][10]={0};
doubleB[10]={0},C[10]={0};
inti,j;
for（i=1;i<=10;i++）
for（j=1;j<=10;j++）
{
if（i!
=j）
A[i-1][j-1]=sin（0.5*i+0.2*j）;
else
A[i-1][j-1]=1.52*cos（i+1.2*j）;
}
//对实矩阵A进行拟上三角化
nishangsanjiaohua（A）;
zhengli（A）;
cout<<"矩阵A经过拟上三角化所得的矩阵A（n-1）:
"<for（i=0;i{
for（j=0;j{
cout<:
scientific）<}
cout<<'\n'<}
//拟上三角化后进行的QR分解
QRfenjie（A,Q,R）;
zhengli（R）;
cout<<"矩阵A（n-1）三角化得到的Q矩阵:
"<for（i=0;i{
for（j=0;j{
cout<:
scientific）<}
cout<<'\n'<}
cout<<"矩阵A（n-1）三角化得到的R矩阵:
"<for（i=0;i{
for（j=0;j{
cout<:
scientific）<}
cout<<'\n'<}
//求解A矩阵的全部特征值
subuQR（A,B,C）;
zhengli（A）;
cout<<"矩阵A（n-1）双步位移QR迭代后RQ阵:
"<for（i=0;i{
for（j=0;j{
cout<:
scientific）<}
cout<<'\n'<}
cout<<"矩阵A（n-1）双步位移QR迭代后求出的所有特征值:
"<for（i=0;i{
cout<:
scientific）<cout<<'\n'<}
for（i=1;i<=10;i++）
for（j=1;j<=10;j++）
{
if（i!
=j）
A[i-1][j-1]=sin（0.5*i+0.2*j）;
else
A[i-1][j-1]=1.52*cos（i+1.2*j）;
}
//A相应于实特征值的特征向量
cout<<"矩阵A（n-1）双步位移QR迭代后求出的所有实特征值所对应的特征向量:
"<for（i=0;i{
if（C[i]==0）
{
cout<<"λ["<tezhengxl（A,B[i]）;
}
}
return1;
}
voidzhengli（double（*A）[10]）
{
inti,j;
for（i=0;ifor（j=0;j{
if（fabs（A[i][j]）<=E）
A[i][j]=0;
}
}
voidnishangsanjiaohua（double（*A）[10]）
{
doubled,c,h,sum,t;
doubleu[N],p[N],q[N],w[N];
intr,i,j,k;
for（r=0;r<=N-3;r++）
{
k=0;
//
（1）做判断
for（i=r+2;i{
if（A[i][r]==0）
k++;
}
if（k!
=N-r-2）
{
//
（2）计算
sum=0;
for（i=r+1;i{
sum=sum+A[i][r]*A[i][r];
}
d=sqrt（sum）;
if（A[r+1][r]==0）
c=d;
else
c=（-1）*fabs（A[r+1][r]）/A[r+1][r]*d;
h=c*c-c*A[r+1][r];
//（3）给u赋值
for（i=0;i{
if（i<=r）
u[i]=0;
else
{
if（i==r+1）
u[i]=A[i][r]-c;
else
u[i]=A[i][r];
}
}
//（4）计算
for（i=0;i{
u[i]=u[i]/h;
}
for（i=0;i{
sum=0;
for（j=r+1;j{
sum=sum+A[j][i]*u[j];
}
p[i]=sum;
sum=0;
for（j=r+1;j{
sum=sum+A[i][j]*u[j];
}
q[i]=sum;
}
sum=0;
for（i=r+1;i{
sum=sum+p[i]*u[i];
}
t=sum;
for（i=0;i{
w[i]=q[i]-t*u[i]*h;
}
for（i=0;ifor（j=0;j{
A[i][j]=A[i][j]-w[i]*u[j]*h-u[i]*h*p[j];
}
}
}
}
voidQRfenjie（double（*A）[10],double（*Q）[N],double（*R）[N]）
{
inti,j,k,m;
doubled,c,h;
doubleu[N],w[N],p[N];
for（i=0;i{
for（j=0;j{
if（i==j）Q[i][j]=1;
elseQ[i][j]=0;
}
}
for（i=0;i{
for（j=0;jR[i][j]=A[i][j];
}
for（i=0;i{
for（j=i+1;jif（R[j][i]<=E）m=m+1;
if（m==（N-1-i））continue;
else
{
for（j=i,d=0;jd=d+R[j][i]*R[j][i];
d=sqrt（d）;
c=-1*fabs（R[i][i]）/R[i][i]*d;
h=c*c-c*R[i][i];
for（j=i+1;ju[j]=R[j][i];
for（j=0;ju[j]=0;
u[i]=R[i][i]-c;
for（j=0;j{
for（k=0,w[j]=0;kw[j]=Q[j][k]*u[k]+w[j];
}
for（j=0;j{
for（k=0;kQ[j][k]=Q[j][k]-w[j]*u[k]/h;
}
for（j=0;j{
for（k=i,p[j]=0;kp[j]=R[k][j]*u[k]+p[j];
p[j]=p[j]/h;
}
for（j=0;j{
for（k=0;kR[j][k]=R[j][k]-u[j]*p[k];
}
}
}
}//矩阵的QR分解
doublekaifang（doubleb,doublec）
{
doublem;
m=b*b-4*c;
returnm;
}
//使用双步位移QR法求实矩阵A的全部特征值
voidsubuQR（double（*A）[10],double*RR,double*II）
{
intm=N-1,BU=3,i,j;
intL=1;
intk=0;
doubles,t,x;
doubleM[N][N],B[N][N];
intf=0;
doubled,c,h;
doubleu[N],w[N],p[N];
doubleQ[N][N],R[N][N];
while（BU!
=11）
{
//编程精妙之处*****************************
if（BU==3）
{
if（fabs（A[m][m-1]）<=E）
{
RR[m]=A[m][m];II[m]=0;
m=m-1;
BU=4;
}
else
BU=5;
}
if（BU==4）
{
if（m==0）
{
RR[m]=A[m][m];
II[m]=0;
BU=11;
}
else
{
BU=3;
}
}
if（BU==5）
{
s=A[m-1][m-1]+A[m][m];
t=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m];
x=kaifang（s,t）;
if（x>=0）
{
x=sqrt（x）;
RR[m]=（s-x）/2;RR[m-1]=（s+x）/2;
II[m]=0;II[m-1]=0;
}
else
{
x=sqrt（-x）;
RR[m]=s/2;RR[m-1]=s/2;
II[m]=x/2;II[m-1]=-x/2;
}
BU=6;
}
if（BU==6）
{
if（m==1）
{
BU=11;
}
else
BU=7;
}
if（BU==7）
{
if（fabs（A[m-1][m-2]）<=E）
{
m=m-2;BU=4;
}
else
BU=8;
}
if（BU==8）
{
if（L==MAX）
BU=11;
else
BU=9;
}
if（BU==9）
{
for（i=0;i<=m;i++）
{
for（j=0;j<=m;j++）
{
if（i==j）
{
for（k=0,M[i][j]=0;k<=m;k++）
M[i][j]=A[i][k]*A[k][j]+M[i][j];
M[i][j]=M[i][j]-s*A[i][j]+t;
}
else
{
for（k=0,M[i][j]=0;k<=m;k++）
M[i][j]=A[i][k]*A[k][j]+M[i][j];
M[i][j]=M[i][j]-s*A[i][j];
}
}
}
//以下是M的QR分解
for（i=0;i<=m;i++）
{for（j=0;j<=m;j++）{if（i==j）Q[i][j]=1;elseQ[i][j]=0;}}
for（i=0;i<=m;i++）
{for（j=0;j<=m;j++）R[i][j]=M[i][j];}
for（i=0;i{
for（j=i+1;j<=m;j++）if（R[j][i]<=E）f=f+1;
if（f==（m-i））continue;
for（j=i,d=0;j<=m;j++）d=d+R[j][i]*R[j][i];
d=sqrt（d）;
c=-1*fabs（R[i][i]）/R[i][i]*d;
h=c*c-c*R[i][i];
for（j=i+1;j<=m;j++）u[j]=R[j][i];
for（j=0;ju[i]=R[i][i]-c;
for（j=0;j<=m;j++）
{for（k=0,w[j]=0;k<=m;k++）w[j]=Q[j][k]*u[k]+w[j];}
for（j=0;j<=m;j++）
{for（k=0;k<=m;k++）Q[j][k]=Q[j][k]-w[j]*u[k]/h;}
for（j=0;j<=m;j++）
{for（k=i,p[j]=0;k<=m;k++）p[j]=R[k][j]*u[k]+p[j];p[j]=p[j]/h;}
for（j=0;j<=m;j++）
{
for（k=0;k<=m;k++）R[j][k]=R[j][k]-u[j]*p[k];
}
}
for（j=0;j<=m;j++）
{
for（k=0;k<=m;k++）M[j][k]=Q[j][k];
}
for（i=0;i<=m;i++）
{
for（j=0;j<=m;j++）
{for（k=0,B[i][j]=0;k<=m;k++）B[i][j]=M[k][i]*A[k][j]+B[i][j];}
}
for（i=0;i<=m;i++）
{
for（j=0;j<=m;j++）
{for（k=0,A[i][j]=0;k<=m;k++）A[i][j]=B[i][k]*M[k][j]+A[i][j];}
}
BU=10;
}
if（BU==10）
{
L=L+1;
BU=3;
}
}
}
voidtezhengxl（double（*a）[N],doubleT）
{
voidqstzxl（double（*a）[N]）;
doubleM[N][N];
inti,j;
for（i=0;ifor（j=0;j{
if（i==j）
M[i][j]=a[i][j]-T;
else
M[i][j]=a[i][j];
}
qstzxl（M）;
}
//巧妙的使用老师上课要求上机调试练习的列主元高斯消去法求解实特征值对应的特征函数
voidqstzxl（double（*a）[N]）
{
doubleb[N]={0};
doubleH[N][N]={0},l[N]={0};
doubleX[N];
doubleB;
doublesum;
inti,j,m,k,z;
for（k=0;k{
for（j=k;j{
l[j]=a[k][j];
}
z=k;
for（m=k;m{
if（fabs（a[z][k]）z=m;
}
for（j=k;j{
a[k][j]=a[z][j];
a[z][j]=l[j];
}
B=b[k];
b[k]=b[z];
b[z]=B;
for（i=k+1;i{
H[i][k]=a[i][k]/a[k][k];
for（j=k+1;ja[i][j]=a[i][j]-H[i][k]*a[k][j];
b[i]=b[i]-H[i][k]*b[k];
}
}
X[N-1]=1;
for（k=N-2;k>=0;k--）
{
sum=0;
for（j=k+1;j{
sum=sum+a[k][j]*X[j];
}
X[k]=（b[k]-sum）/a[k][k];
}
cout<<"（";
for（i=0;icout<cout<<"）";
cout<<'\n'<}
3、
程序运行结果显示
四、现象讨论和分析
由于我们知道矩阵A和拟上三角矩阵A（n-1）是相似矩阵，具有相同的特征向量和特征值。
在程序调试的过程当中发现，在求实值特征值对应的特征向量时，分别使用矩阵A和拟上三角矩阵A（n-1）得到特征向量是有较大区别的如图3所示。
由于拟上三角矩阵左下角0较多，而且在减去特征值对角矩阵时，会把误差积累增加。
因此建议使用原矩阵A来计算实特征值对应的特征向量。
升级会员 