北航数值分析第二次大作业.docx
《北航数值分析第二次大作业.docx》由会员分享，可在线阅读，更多相关《北航数值分析第二次大作业.docx（26页珍藏版）》请在冰豆网上搜索。
北航数值分析第二次大作业
数值分析第二次大作业
1、算法的设计方案
首先构造矩阵A，利用Householder矩阵对矩阵A作相似变换，把A化为拟上三角矩阵A（n-1），其算法如书上第61页所示。
使用QR分解法对矩阵A（n-1）进行QR分解，其算法如书上第59页所示，进而可以求出所得矩阵的Q、R、RQ阵。
然后对A（n-1）进行带双步位移的QR分解求解矩阵的全部特征值。
这里采用这么几步进行：
第一步：
判断是否
，若不是，则进入第四步。
若是，则
为特征值，m=m-1，若m=1，进入第二步，若m=2进入第三步，否则进入第四步。
第二步：
m=1，则
为特征值，转向结束步。
第三步：
m=2，则可以求出A的两个特征值
和
，转向结束步。
第四步：
判断是否
，若不是，进入第五步。
若是，则得到A的两个特征值
和
，令m=m-2，若m=1，进入第二步，若m=2进入第三步，否则进入第一步。
第五步：
判断是否达到循环上限，若达到上限，则结束，否则进入第六步。
第六步：
对A进行双步移位QR分解，这里的算法如书上第64页所示，分解之后转向计数步。
计数步：
对循环次数进行计数，并转向第一步。
结束步：
显示所求得的特征值。
最后对实特征值利用列主元高斯消元法求解其对应的特征向量，其算法见书第17页。
2、对程序的几点说明
1、对A求拟上三角化矩阵A（n-1），这里采用一个单独的函数NSSjiao（）来求出。
2、对矩阵A（n-1）进行QR分解采用一个单独的函数QRfenjie（）求出。
3、对矩阵A（n-1）进行带双步位移的QR分解采用一个单独的函数SBQR（）求出。
4、对实特征值利用列主元高斯消元法求解其对应的特征向量在这里采用一个单独的函数vector（）来实现。
5、所有的结果都用C语言文件指令使之显示在一个txt文件中，方便结果的复制与粘贴。
A（n-1）输出到NSSjiao.txt中，A（n-1）QR分解所得矩阵的Q、R、RQ阵输出到q.txt中，A的特征值输出到tezhenzhi.txt中，A实特征值对应的特征向量输出到vector.txt中。
三、源程序如下：
#include"stdio.h"
#include"math.h"
#defineN10+1
#definen10
#defineL1000
doublea[N][N];
doubleRe[N]={0};
doubleIm[N]={0};
intsgn（doublea）;//符号函数
voidNSSjiao（）;//拟上三角化函数
voidQRfenjie（）;//QR分解函数
voidSBQR（）;//双步QR分解函数
voidvector（）;//求特征向量的函数
voidmain（）
{
inti,j;
for（i=1;i<=n;i++）
{for（j=1;j<=n;j++）
{if（i!
=j）
a[i][j]=sin（0.5*i+0.2*j）;
else
a[i][j]=1.5*cos（i+1.2*j）;
}
}//赋值，矩阵初始化
NSSjiao（）;
QRfenjie（）;
SBQR（）;
vector（）;
}
intsgn（doublea）
{if（a>=0）
return
（1）;
elsereturn（-1）;
}
voidNSSjiao（）
{intr,i,j;
intcnt=0;
doubledr=0;
doublecr=0;
doublehr=0;
doubletr=0;
doubleur[N]={0};
doublepr[N]={0};
doubleqr[N]={0};
doublewr[N]={0};
for（r=1;r<=n-2;r++）
{for（i=r+2;i<=n;i++）
{
if（a[i][r]==0）
cnt++;
}
if（cnt!
=n-r-1）
{for（i=r+1;i<=n;i++）
dr+=a[i][r]*a[i][r];
dr=sqrt（dr）;
cr=-sgn（a[r+1][r]）*dr;
hr=cr*cr-cr*a[r+1][r];
for（i=1;i<=r;i++）
ur[i]=0;
for（i=r+2;i<=n;i++）
ur[i]=a[i][r];
ur[r+1]=a[r+1][r]-cr;
for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
pr[i]+=a[j][i]*ur[j]/hr;
for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
qr[i]+=a[i][j]*ur[j]/hr;
for（i=1;i<=n;i++）
tr+=pr[i]*ur[i]/hr;
for（i=1;i<=n;i++）
wr[i]=qr[i]-tr*ur[i];
for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
a[i][j]-=（wr[i]*ur[j]+ur[i]*pr[j]）;
dr=0;
cr=0;
hr=0;
tr=0;
for（i=1;i<=n;i++）
{
pr[i]=0;
qr[i]=0;
}
}
}
FILE*fp;
fp=fopen（"NSSjiao.txt","a"）;
fprintf（fp,"A矩阵拟上三角化后为：
\n"）;
for（i=1;i<=n;i++）
{
for（j=1;j<=n;j++）
fprintf（fp,"%13.11e\t",a[i][j]）;
fprintf（fp,"\n"）;
}
fclose（fp）;
}
voidQRfenjie（）
{intr,i,j,t;
intcnt=0;
doubledr=0;
doublecr=0;
doublehr=0;
doubleur[N]={0};
doublewr[N]={0};
doublepr[N]={0};
doubleQ[N][N]={0};
doubleRQ[N][N]={0};
for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
{if（i==j）
Q[i][j]=1;//Q阵初始化
}
for（r=1;r<=n-1;r++）
{for（i=r+1;i<=n;i++）
{if（a[i][r]==0）
cnt++;
}
if（cnt!
=n-r）
{for（i=r;i<=n;i++）
dr+=a[i][r]*a[i][r];
dr=sqrt（dr）;
cr=-sgn（a[r][r]）*dr;
hr=cr*cr-cr*a[r][r];
for（i=1;iur[i]=0;
for（i=r+1;i<=n;i++）
ur[i]=a[i][r];
ur[r]=a[r][r]-cr;
for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
wr[i]+=Q[i][j]*ur[j];
for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
Q[i][j]-=wr[i]*ur[j]/hr;
for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
pr[i]+=a[j][i]*ur[j]/hr;
for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
a[i][j]-=ur[i]*pr[j];
dr=0;
cr=0;
hr=0;
for（i=1;i<=n;i++）
{
wr[i]=0;
pr[i]=0;
}
}
}
FILE*fp;
fp=fopen（"q.txt","a"）;
fprintf（fp,"A（n-1）矩阵QR分解后Q阵为:
\n"）;
for（i=1;i<=n;i++）
{
for（j=1;j<=n;j++）
fprintf（fp,"%13.11e\t",Q[i][j]）;
fprintf（fp,"\n"）;
}
fclose（fp）;
fp=fopen（"Q.txt","a"）;
fprintf（fp,"A（n-1）矩阵QR分解后R阵为：
\n"）;
for（i=1;i<=n;i++）
{
for（j=1;j<=n;j++）
fprintf（fp,"%13.11e\t",a[i][j]）;
fprintf（fp,"\n"）;
}
fclose（fp）;
fp=fopen（"Q.txt","a"）;
fprintf（fp,"A（n-1）矩阵QR分解后RQ阵为：
\n"）;
for（i=1;i<=n;i++）
{for（j=1;j<=n;j++）
{for（t=1;t<=n;t++）
RQ[i][j]+=a[i][t]*Q[t][j];
fprintf（fp,"%13.11e\t",RQ[i][j]）;
}
fprintf（fp,"\n"）;
}
fclose（fp）;
}
voidSBQR（）
{inti,j;
intm,k;
intr;
intcnt;
intcount;
doubles=0,t=0;
doubledelta=0;
doubledr=0;
doublecr=0;
doublehr=0;
doubletr=0;
doubleur[N]={0};
doublewr[N]={0};
doublepr[N]={0};
doublevr[N]={0};
doubleqr[N]={0};
doubleI[N][N]={0};
doubleaa[N][N]={0};
doublec[N][N]={0};
doubleb[N][N]={0};
doubleM[N][N]={0};
for（i=1;i<=n;i++）
{for（j=1;j<=n;j++）
{if（i!
=j）
a[i][j]=sin（0.5*i+0.2*j）;
else
a[i][j]=1.5*cos（i+1.2*j）;
}
}//赋值，矩阵初始化
NSSjiao（）;
m=n;
count=1;
gotoloop1;
loop1:
if（fabs（a[m][m-1]）<=1e-12）//降一阶的判断
{
Re[m]=a[m][m];
Im[m]=0;
m--;
if（m==1）
gotoloop2;
if（m==2）
gotoloop3;
else
gotoloop4;
}
else
gotoloop4;
loop2:
//m=1
Re[m]=a[1][1];
Im[m]=0;
gotoloopend;
loop3:
//m=2
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];
delta=s*s-4*t;
if（delta<=0）
{
Re[m]=s/2;
Im[m]=sqrt（fabs（delta））/2;
m--;
Re[m]=s/2;
Im[m]=sqrt（fabs（delta））/2;
m--;
}
else
{Re[m]=s/2+sqrt（fabs（delta））/2;
Im[m]=0;
m--;
Re[m]=-s/2-sqrt（fabs（delta））/2;
Im[m]=0;
m--;
}
gotoloopend;
loop4:
if（fabs（a[m-1][m-2]）<=1e-12）//降二阶的判断
{
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];
delta=s*s-4*t;
if（delta<=0）
{
Re[m]=s/2;
Im[m]=sqrt（fabs（delta））/2;
m--;
Re[m]=s/2;
Im[m]=-sqrt（fabs（delta））/2;
m--;
}
else
{Re[m]=s/2+sqrt（fabs（delta））/2;
Im[m]=0;
m--;
Re[m]=s/2-sqrt（fabs（delta））/2;
Im[m]=0;
m--;
}
if（m==1）
gotoloop2;
elseif（m==2）
gotoloop3;
elsegotoloop1;
}
elsegotoloop5;
loop5:
if（count==L）
gotoloopend;
elsegotoloop6;
loop6:
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];
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）
{if（i==j）
I[i][j]=1;
}
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）
{for（k=1;k<=m;k++）
aa[i][j]+=a[i][k]*a[k][j];
}
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）
{M[i][j]=aa[i][j]-s*a[i][j]+t*I[i][j];
b[i][j]=M[i][j];
c[i][j]=a[i][j];
}
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）
aa[i][j]=0;
for（r=1;r<=m-1;r++）
{for（i=r+1;i<=m;i++）
{if（b[i][r]==0）
cnt++;
}
if（cnt!
=m-r）
{for（i=r;i<=m;i++）
dr+=b[i][r]*b[i][r];
dr=sqrt（dr）;
cr=-sgn（b[r][r]）*dr;
hr=cr*cr-cr*b[r][r];
for（i=1;iur[i]=0;
for（i=r+1;i<=m;i++）
ur[i]=b[i][r];
ur[r]=b[r][r]-cr;
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）
vr[i]+=b[j][i]*ur[j]/hr;
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）
b[i][j]-=ur[i]*vr[j];
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）
pr[i]+=c[j][i]*ur[j]/hr;
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）
qr[i]+=c[i][j]*ur[j]/hr;
for（i=1;i<=m;i++）
tr+=pr[i]*ur[i]/hr;
for（i=1;i<=m;i++）
wr[i]=qr[i]-tr*ur[i];
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）
c[i][j]-=（wr[i]*ur[j]+ur[i]*pr[j]）;
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）
a[i][j]=c[i][j];
dr=0;
cr=0;
hr=0;
tr=0;
for（i=1;i<=m;i++）
{
pr[i]=0;
qr[i]=0;
vr[i]=0;
}
}
}
gotoloopcount;
loopcount:
count++;
gotoloop1;
loopend:
FILE*fp;
fp=fopen（"tezhenzhi.txt","a"）;
fprintf（fp,"矩阵A的特征值为:
\n"）;
for（i=1;i<=10;i++）
{
if（Im[i]>=0）
fprintf（fp,"%13.11e+%13.11e*i\n",Re[i],Im[i]）;
else
fprintf（fp,"%13.11e%13.11e*i\n",Re[i],Im[i]）;
}
fclose（fp）;
}
voidvector（）//列主元高斯消元法求解特征向量
{intcount;
intk,i,j,ik;
doublem=0;
doublex[N]={0};
doubletemp=0;
doublesum=0;
for（count=1;count<=n;count++）
{if（Im[count]==0）
{
for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
{
if（i!
=j）
a[i][j]=sin（0.5*i+0.2*j）;
else
a[i][j]=1.5*cos（i+1.2*j）-Re[count];
}
for（k=1;k<=n;k++）
{
ik=k;
for（i=k;i<=n;i++）
{
if（fabs（a[i][k]）>fabs（a[ik][k]））
{
ik=i;
}
}
for（j=k;j<=n;j++）
{
temp=a[k][j];
a[k][j]=a[ik][j];
a[ik][j]=temp;
}
for（i=k+1;i<=n;i++）
{
m=a[i][k]/a[k][k];
for（j=k+1;j<=10;j++）
a[i][j]=a[i][j]-m*a[k][j];
}
}
x[n]=1;
for（k=n-1;k>=1;k--）
{
temp=0;
for（j=k+1;j<=n;j++）
{
temp+=a[k][j]*x[j];
x[k]=-temp/a[k][k];
}
}
FILE*fp;
fp=fopen（"vector.txt","a"）;
fprintf（fp,"属于特征值%13.11e的特征向量为\n",Re[count]）;
for（i=1;i<=n;i++）
fprintf（fp,"%13.11e\t",x[i]）;
fprintf（fp,"\n",0）;
fclose（fp）;
}
}
}
3、运行结果如下
对A进行拟上三角化得到的矩阵A（n-1）:
-8.82751675883e-001-9.93313649183e-002-1.10334928599e+000-7.60044358564e-0011.54910107991e-001-1.94659186287e+000-8.78243638293e-002-9.25588938718e-0016.03259944053e-0011.51886095647e-001
-2.34787836242e+0002.37237010494e+0001.81929082221e+0003.23780410155e-0012.20579844032e-0012.10269266255e+0001.81613808610e-0011.27883908999e+000-6.38057812440e-001-4.15407560380e-001
0.00000000000e+0001.72827459997e+000-1.17146764279e+000-1.24383926270e+000-6.39975834174e-001-2.00283307904e+0002.92494720612e-001-6.41283006839e-0019.78399762128e-0022.55776357416e-001
0.00000000000e+0000.00000000000e+000-1.29166953413e+000-1.11160351340e+0001.17134682410e+000-1.30735603002e+0001.80369917775e-001-4.24638535837e-0017.98895523930e-0021.60881992807e-001
0.00000000000e+0000.00000000000e+000-8.29733316381
升级会员 