北航数值分析实习第二题.docx
《北航数值分析实习第二题.docx》由会员分享，可在线阅读，更多相关《北航数值分析实习第二题.docx（17页珍藏版）》请在冰豆网上搜索。
数值分析计算实习报告
第二题
所在班级
A1班
学生姓名
学生学号
2015年11月
北航数值分析计算实习报告 第16页
1算法设计方案
1.1矩阵的输入
A矩阵是一个10乘10的矩阵，同时并没有对称性，所以不能压缩，直接存储即可。
1.2对矩阵进行上三角化
A矩阵的上三角化采用常规方法即可，得到。
1.3QR分解法
对上三角化后的矩阵进行QR分解可以减少计算量。
RQ这个矩阵是与原矩阵A相似的矩阵，与A具有相同的特征值和特征向量。
1.4求解特征值
采用带双步位移的QR分解法对求解特征值即可。
1.5求解实特征值的特征向量
求解特征向量，本质上就是求解方程
的解，其中I是单位向量。
直接采用列主元素Gauss消去法求解即可。
但是，值得注意的是，对应于一个特征值的特征向量是无穷多的，也就是经过Gauss消去法后，最后一行全是0。
因此，求解时候，需要给特征向量某个元素赋值，本题，将最后一个元素统一赋值为1，即可求出对应于某个特征值的特征向量的基础解系，而全部特征向量为，k取不等于0的任意实数。
2C++程序
#include
#include
#include
#include
#definen11
#defineerr1e-12
#defineL2500
voidcaculateA（）;//计算矩阵A的系数
voidhessenberg（）;//将A拟上三角化
voidQR（）;//对矩阵进行QR分解
voidQRshuangbu（）;//对矩阵进行带双步位移的QR分解
voidgauss（）;//列主元的高斯消元法求解特征向量
inti,j,s,p,k,ik,nR,nC;
doubleA[n][n],q[n][n],r[n][n],rq[n][n],I[n][n];
doubleP[n],W[n],u[n],Q[n];
doubledr,cr,hr,ar,tr;
doubles1,t,x,tzR[n],tzC[2][n],sum,M[n][n],v[n];
voidmain（）
{for（i=1;i {for（j=1;j {if（i==j）
{I[i][j]=1;q[i][j]=1;}
else{I[i][j]=0;q[i][j]=0;}}}
caculateA（）;
hessenberg（）;
QR（）;
QRshuangbu（）;
gauss（）;
printf（"运行时间为%d\n",clock（））;
}
voidcaculateA（）//计算矩阵A的系数
{
for（i=1;i{for（j=1;j {if（j!
=i）
A[i][j]=sin（0.5*i+0.2*j）;
else
A[i][j]=1.52*cos（i+1.2*j）;}
}
}
voidhessenberg（）//将A拟上三角化
{
for（s=1;s{for（ar=0.0,i=s+2;i ar+=A[i][s]*A[i][s];
if（ar<=err）
continue;
else{
ar+=A[s+1][s]*A[s+1][s];
dr=sqrt（ar）;
if（A[s+1][s]>0）
cr=-dr;
else
cr=dr;
hr=cr*cr-cr*A[s+1][s];
for（i=1;i<=s;i++）
u[i]=0.0;
u[s+1]=A[s+1][s]-cr;
for（i=s+2;i u[i]=A[i][s];
for（j=1;j {for（P[j]=0.0,i=1;i P[j]+=A[i][j]*u[i]/hr;}
for（tr=0.0,i=1;i {tr+=P[i]*u[i]/hr;}
for（i=1;i {for（Q[i]=0.0,j=1;j Q[i]+=A[i][j]*u[j]/hr;}
for（i=1;i {W[i]=Q[i]-tr*u[i];}
for（i=1;i {for（j=1;j A[i][j]-=（W[i]*u[j]+u[i]*P[j]）;}
}
}
for（i=1;i{for（j=1;jif（fabs（A[i][j]） A[i][j]=0.0;}
printf（"数值分析计算实习第二次作业\n"）;
printf（"A1班\n"）;
printf（"\n"）;
printf（"矩阵A经过拟上三角化后得到的矩阵为:
\n"）;
for（i=1;i {
for（j=1;jprintf（"%1.12e,",A[i][j]）;
printf（"\n"）;}
}
voidQR（）//对矩阵进行QR分解
{
doubleu[n],w[n],F[n];
for（s=1;s for（p=1;p r[s][p]=A[s][p];
for（s=1;s {
for（dr=0.0,i=s;i dr+=r[i][s]*r[i][s];
dr=sqrt（dr）;
if（dr<=err）
continue;
else
{if（A[s][s]>0）cr=-dr;
elsecr=dr;
hr=cr*cr-cr*r[s][s];
for（i=1;i u[i]=0;
u[s]=r[s][s]-cr;
for（i=s+1;i u[i]=r[i][s];
for（i=1;i {for（F[i]=0.0,j=1;j F[i]+=r[j][i]*u[j]/hr;}
for（i=1;i {for（j=1;j r[i][j]=r[i][j]-u[i]*F[j];}
for（i=1;i {for（w[i]=0.0,j=1;j w[i]+=q[i][j]*u[j];}
for（i=1;i {for（j=1;j q[i][j]-=w[i]*u[j]/hr;}
}
}
for（i=1;i {for（j=1;j {for（rq[i][j]=0.0,s=1;s rq[i][j]+=r[i][s]*q[s][j];}}
printf（"\n"）;
printf（"A经过QR分解后得到的正交矩阵Q为：
\n"）;
for（i=1;i{for（j=1;jif（fabs（q[i][j]）q[i][j]=0.0;}
for（i=1;i {for（j=1;jprintf（"%1.12e,",q[i][j]）;
printf（"\n"）;}
printf（"\n"）;
printf（"A经过QR分解后得到的上三角矩阵R为：
\n"）;
for（i=1;i{for（j=1;jif（fabs（r[i][j]）r[i][j]=0.0;}
for（i=1;i {for（j=1;jprintf（"%1.12e,",r[i][j]）;
printf（"\n"）;}
printf（"\n"）;
printf（"正交矩阵Q乘以上三角矩阵R得到的RQ为：
\n"）;
for（i=1;i{for（j=1;jif（fabs（rq[i][j]）rq[i][j]=0.0;}
for（i=1;i {for（j=1;jprintf（"%1.12e,",rq[i][j]）;
printf（"\n"）;}
}
voidQRshuangbu（）//对矩阵进行带双步位移的QR分解
{
intK=1,m=10;
nR=0,nC=0;
loop3:
if（fabs（A[m][m-1]）<=err）
{nR++;
tzR[nR]=A[m][m];
m--;
gotoloop4;
}
elsegotoloop5;
loop4:
if（m==1）
{
nR++;
tzR[nR]=A[1][1];
gotoloop9;
}
elseif（m==2）
{
s1=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=s1*s1-4*t;
if（x>=0）
{
nR++;
tzR[nR]=（s1+sqrt（x））/2;
nR++;
tzR[nR]=（s1-sqrt（x））/2;
}
else
{
nC++;
tzC[0][nC]=s1/2;
tzC[1][nC]=sqrt（-x）/2;
nC++;
tzC[0][nC]=s1/2;
tzC[1][nC]=-sqrt（-x）/2;
}
gotoloop9;
}
elsegotoloop3;
loop5:
{
if（fabs（A[m-1][m-2]）<=err）
{
s1=A[m
升级会员 