北航数值分析大作业2.docx
《北航数值分析大作业2.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业2.docx（23页珍藏版）》请在冰豆网上搜索。
北航数值分析大作业2
一、算法设计方案
设计思想：
用带双步位移的QR分解法求矩阵A（10*10）的全部特征值。
在计算出A的基础上，先利用Householder矩阵对矩阵A作相似变换，把A化为拟上三角矩阵A（n-1），然后进行带双步位移的QR分解（其中Mk的QR分解可调用子程序），通过调用一元二次方程的根解二阶块矩阵的特征值，最后计算出A（n-1）的特征值，即为A的特征值，然后对实特征值利用列主元高斯消元法求解其对应的特征向量。
具体算法如下：
（给定精度Epsilon=le-20，最大迭代次数L=1000,n=10,N=11）
1.计算矩阵A
对i,j=1到10执行
（1）ifi=j,a[i][j]=1.5*cos*（i+1.2j）;elsea[i][j]=sin*（0.5*i+0.2j）;
（2）prinfta[i][j].
2.对矩阵A拟上三角化A（n-1）
A[i][j]为矩阵A的元素,i,j=1,…,n。
设置标志量flag=0，flag=0表示A[i][r]全为零，flag=1表示A[i][r]不全为零。
对r=1到n-2执行
（1）if|air|>Epsilon,flag=1,转
（2）；
（2）计算
dr=sqrt（A[i][r]*A[i][r]）
cr=--sgn（A[r+1][r]）
hr=cr*（cr-A[r+1][r]）;
（3）对i=r+1到n执行ur[i]=A[i][r];ur[r+1]=ur[r+1]-cr；
（4）计算
pr=A’*ur/hr
qr=A*ur/hr
tr=pr’*ur/hr
wr=qr-tr*ur
A=A-wr*ur-ur*pr’
3.对A（n-1）做双步位移QR分解
（1）记A1=A（n-1）=A[i][j]，令k=0,m=n;
（2）if|A[m][m-1]<=Epsilon|,则得到A的一个特征值Lamda[m][0]=A[m][m]，置m=m-1,转（3）;否则转（4）;
（3）ifm=1,则得到A的一个特征值Lamda[m][0]=A[1][1]，转（10）;ifm=0，则直接转（10）；ifm>1,则转
（2）;
（4）求二阶子阵的两个特征值，即计算二次方程的两个根，
x*x-s1[0]*x+s2[0]=0（*）
其中s1[0]=A[m-1][m-1]+A[m][m];s2[0]=A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1]
（该过程可用求解二次方程子程序实现，注意讨论根的情况）；
（5）ifm=2,则得到A的两个特征值，转（10）；否则转（7）；
（6）if|A[m-1][m-2]）|<=Epsilon,则得到A的两个特征根，置m=m-2，转（3）;否则转（7）；
（7）ifk=L,printf“FailtogetalltheeigenvaluesofmatrixA!
”；否则转（8）；
（8）对i,j=1,…,m，计算
s=A[m-1][m-1]+A[m][m];
t=A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1];
Mk[i][j]=Mk[i][j]-s*A[i][j]+t*I[i][j]
Mk=Qk*Rk（对Mk作QR分解，可调用QR分解子程序）
Ak+1=Qk’*Ak*Qk
（9）置k=k+1，转
（2）；
（10）printf“GetalltheeigenvaluesofmatrixA!
Congratulations!
”，停止计算。
Mk的QR分解算法：
记B[][N]=Mk[][N],C[][N]=A[][N]，置flag=0（B[i][r]=0）。
对r=1到m-1执行
（1）if|B[i][r]|>=Epsilon,flag=1，转
（2）;
（2）计算
dr=sqrt（A[i][r]*A[i][r]）
cr=--sgn（A[r+1][r]）
hr=cr*（cr-A[r+1][r]）;
（3）对i=r+1到m执行ur[i]=A[i][r];ur[r+1]=ur[r+1]-cr；
（4）计算
vr=Br’*ur/hr
Br+1=Br-ur*vr’
pr=Crur/hr
qr=Crr/hr
tr=pr’*ur/hr
wr=qr-tr*ur
Cr+1=Cr-wr*ur-ur*pr’
（5）continue.
求二元一次方程算法：
置根Lamda[][3]={0}。
对方程（*），令t1=s1[0];t2=s2[0];
计算
delt=t1*t1-4*t2;
sqdelt=sqrt（fabs（delt））;
（2）ifdelt>=0.0，则s1[0]=（t1+sqdelt）/2.0，s2[0]=（t1-sqdelt）/2.0;否则转（3）;
（3）置Lamda[m][2]=1;Lamda[m-1][2]=1;s1[0]=t1/2.0;s1[1]=sqdelt/2.0;s2[0]=t1/2.0;s2[1]=-sqdelt/2.0;
4.用列主元Gauss消去法求解实特征根对应的特征向量
（1）ifLamda[m][2]=0，转
（2）;
（2）调用GaussOptimal子程序；
（3）printf“Whentheeigenvalueis%4.12e,thecorrespondingeigenvectoris\n",Lamda[i][0]”
GaussOptimal算法
记a1=aij，令b[N]={0}。
对i=1到10执行A1[i][i]=A1[i][i]-Lamda[p][0]
消元过程：
对于k=1到n-1执行
（1）选行号ik,使得|aik,k|=max|aik|;
（2）交换akj与aik,j（j=1到n）；
（3）对于i=k+1到n计算
mik=A1[i][k]/A1[k][k]
A1[i][j]=A1[i][j]-mik*A1[k][j]（j=k+1到n）
回代过程
（1）iffabs（A1[n][n]）<=1e-9则iffabs（A1[n-1][n-1]）<=1e-9，printf（"两重根\n"），xm=（0,0,…,0,1）,xm-1=（*,*,…,*,1,0）;否则转（3）；
（2）对k=n-2到1执行
1）置temp=0.0;
2）对j=k+1到n执行：
temp=temp+A1[k][j]*x[p][j];
3）x[p-1][k]=-temp/A1[k][k];
（3）printf（"\n单根"）;x[p][n]=1.0;对k=n-1到1执行
1）temp=0.0;
2）对j=k+1到n执行temp=temp+A1[k][j]*x[p][j];
3）x[p][k]=-temp/A1[k][k];
（4）特征向量归一化.
二、源程序
#include/*定义头文件*/
#include
#include
#defineEpsilon1e-20/*定义精度*/
#definen10/*定义矩阵和向量的长度*/
#defineN11/*定义数组用，11时和日常表达数组相符合*/
#defineL1000/*定义最大迭代次数*/
voidHessenberg（doubleA[][N],doubleP[][N]）;/*拟上三角化程序*/
voidQRreduce（doubleC[][N],doubleB[][N],intm）;/*QR分解法*/
intDoublestep（doubleA[][N],doubleLamda[][3],FILE*fp）;/*双步位移的QR分解法*/
voidroot（doubles1[],doubles2[],doubleLamda[][3],intm）;/*求二元一次方程的根*/
voidGaussOptimal（doublea1[][N],doublex[][N],doubleLamda[][3],intp）;/*高斯消元法*/
voidEigenvector（doubleA[][N],doublex[][N],doubleLamda[][3],FILE*fp）;/*求解特征向量对应的特征根*/
intmain（）
{
inti=0,j=0,flag=0;
doublea[N][N]={0.0},P[N][N]={0.0},Lamda[N][3]={0.0},x[N][N]={0.0},a1[N][N]={0.0};
FILE*fp;
fp=fopen（"QRdata.txt","wt"）;/*以只写方式打开QR-data文本文档*/
printf（"*****************************************************\n"）;
printf（"用带双步位移的QR分解方法求解矩阵（10*10）的全部特征值\n"）;/*输出标题*/
printf（"*****************************************************\n"）;
/*计算矩阵A*/
for（i=1;i<=n;i++）
{
for（j=1;j<=n;j++）
{
if（i==j）
{
a[i][j]=1.5*cos（i+1.2*j）;
P[i][j]=1.0;
}
else
a[i][j]=sin（0.5*i+0.2*j）;
}
}
/*输出矩阵A*/
fprintf（fp,"\n原始矩阵A\n"）;
printf（"\n原始矩阵A\n"）;
for（i=1;i<=n;i++）
{
for（j=1;j<=n;j++）
{
fprintf（fp,"%4.12e",a[i][j]）;
printf（"%4.12e",a[i][j]）;
}
fprintf（fp,"\n"）;
printf（"\n"）;
}
printf（"\n********************************************\n"）;
fprintf（fp,"\n********************************************\n"）;
Hessenberg（a,P）;/*拟上三角化*/
printf（"HessenbergisOK!
\n"）;
fprintf（fp,"\n拟上三角化矩阵A（n-1）\n"）;
printf（"\n拟上三角化矩阵A（n-1）\n"）;
for（i=1;i<=n;i++）
{
for（j=1;j<=n;j++）
{
a1[i][j]=a[i][j];
fprintf（fp,"%4.12e",a[i][j]）;
printf（"%4.12e",a[i][j]）;
}
fprintf（fp,"\n"）;
printf（"\n"）;
}
fprintf（fp,"\n拟上三角化变换矩阵P\n"）;
printf（"\n拟上三角化变换矩阵P\n"）;
for（i=1;i<=n;i++）
{
for（j=1;j<=n;j++）
{
printf（"%4.12e",P[i][j]）;
fprintf（fp,"%4.12e",P[i][j]）;
}
fprintf（fp,"\n"）;
printf（"\n"）;
}
printf（"\n********************************************\n"）;
fprintf（fp,"\n********************************************\n"）;
flag=Doublestep（a,Lamda,fp）;/*带双步位移QR分解*/
printf（"\nDoublestepisOK!
\n"）;
printf（"\n矩阵A（n-1）进行QR分解后所得矩阵\n"）;
fprintf（fp,"\n矩阵A（n-1）进行QR分解后所得矩阵\n"）;
for（i=1;i<=N;i++）
{
for（j=1;j<=N;j++）
{
printf（"%4.12e",a[i][j]）;
fprintf（fp,"%4.12e",a[i][j]）;
}
printf（"\n"）;
fprintf（fp,"\n"）;
}
printf（"\n********************************************\n"）;
fprintf（fp,"\n********************************************\n"）;
if（flag）
{
fprintf（fp,"\nTheeigenvaluesofmatrixA\n"）;
printf（"\nTheeigenvaluesofmatrixA\n"）;
for（i=1;i<=n;i++）
{
printf（"%4.12e%4.12e%4.12e\n",Lamda[i][0],Lamda[i][1],Lamda[i][2]）;
fprintf（fp,"%4.12e%4.12e%4.12e\n",Lamda[i][0],Lamda[i][1],Lamda[i][2]）;
}
printf（"\n********************************************\n"）;
fprintf（fp,"\n********************************************\n"）;
}
else
{
printf（"\nFailtogettheeigenvaluesofmatrixA.Pleasecheckit!
\n"）;
return（0）;
}
Eigenvector（a1,x,Lamda,fp）;/*求实特征值对应特征向量*/
printf（"\nEigenvectorisOK!
Getalleigenvectorsoftherealeigenvalues!
\n"）;
fclose（fp）;
return
（1）;
}
/*****************拟上三角化子程序*********************/
voidHessenberg（doubleA[][N],doubleP[][N]）
{
inti,j,r;
intflag=0;/*flag=0全为零，flag=1不全为零*/
doublecr,dr,hr,tr,temp;
doublepr[N]={0.0},qr[N]={0.0},ur[N]={0.0},wr[N]={0.0},p[N]={0.0};
for（r=1;r<=n-2;r++）
{
cr=0.0;dr=0.0;hr=0.0;tr=0.0;temp=0.0;
for（i=1;i<=n;i++）
{
pr[i]=0;qr[i]=0;ur[i]=0;wr[i]=0;p[i]=0;
}
for（i=r+2;i<=n;i++）/*判断第r列是否全为，若有一个不是则flag=1*/
{
if（fabs（A[i][r]）>=Epsilon）
{
flag=1;
break;
}
}
if（flag==1）/*计算dr，cr，hr；给ur赋值*/
{
for（i=r+1;i<=n;i++）
temp=temp+A[i][r]*A[i][r];
dr=sqrt（temp）;
if（A[r+1][r]>Epsilon）
cr=-dr;
else
cr=dr;
hr=cr*（cr-A[r+1][r]）;
for（i=r+1;i<=n;i++）
ur[i]=A[i][r];
ur[r+1]=ur[r+1]-cr;
for（i=1;i<=n;i++）/*计算pr，qr，tr，wr和A*/
{
for（j=1;j<=n;j++）
{
pr[i]=pr[i]+A[j][i]*ur[j];
qr[i]=qr[i]+A[i][j]*ur[j];
p[i]=p[i]+P[i][j]*ur[j];
}
pr[i]=pr[i]/hr;
qr[i]=qr[i]/hr;
p[i]=p[i]/hr;
}
for（i=1;i<=n;i++）
tr=tr+pr[i]*ur[i];
tr=tr/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++）
{
if（（j==r）&&（i>j+1））/*ur前r个元素为零*/
A[i][j]=0.0;
else
A[i][j]=A[i][j]-wr[i]*ur[j]-ur[i]*pr[j];
P[i][j]=P[i][j]-p[i]*ur[j];
}
}
}
}
/*****************QR分解子程序*********************/
voidQRreduce（doubleC[][N],doubleB[][N],intm）
{
inti,j,r;
intflag=0;
doublecr,dr,hr,tr,temp;
doublepr[N],qr[N],ur[N],vr[N],wr[N];
for（r=1;r<=m-1;r++）
{
cr=0.0;dr=0.0;hr=0.0;tr=0.0;temp=0.0;
for（i=1;i<=m;i++）
{
pr[i]=0.0;qr[i]=0.0;ur[i]=0.0;vr[i]=0.0;wr[i]=0.0;
}
for（i=r+1;i<=m;i++）
{
if（fabs（B[i][r]）>=Epsilon）
{
flag=1;
break;
}
}
if（flag==1）
{
for（i=r;i<=m;i++）
temp=temp+B[i][r]*B[i][r];
dr=sqrt（temp）;
if（B[r][r]>Epsilon）
cr=-dr;
else
cr=dr;
hr=cr*（cr-B[r][r]）;
for（i=r;i<=m;i++）
ur[i]=B[i][r];
ur[r]=ur[r]-cr;
for（i=1;i<=m;i++）
{
for（j=1;j<=m;j++）
{
vr[i]=vr[i]+B[j][i]*ur[j];
pr[i]=pr[i]+C[j][i]*ur[j];
qr[i]=qr[i]+C[i][j]*ur[j];
}
vr[i]=vr[i]/hr;
pr[i]=pr[i]/hr;
qr[i]=qr[i]/hr;
}
for（i=1;i<=m;i++）
tr=tr+pr[i]*ur[i];
tr=tr/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]=C[i][j]-wr[i]*ur[j]-ur[i]*pr[j];
B[i][j]=B[i][j]-ur[i]*vr[j];
}
}
}
}
/****************带双步位移QR分解子程序****************/
intDoublestep（doubleA[][N],doubleLamda[][3],FILE*fp）
{
inti,j,k,m,l;
intflag=1;
doubles1[2]={0.0},s2[2]={0.0},s,t;
doubleMk[N][N]={0.0},I[N][N]={0.0};
k=0;m=n;
for（i=1;i<=n;i++）
I[i][i]=1.0;
while
（1）
{
Loop1:
if（fabs（A[m][m-1]）<=Epsilon）
{
Lamda[m][0]=A[m][m];
m=m-1;
Loop2:
if（m==1）
{
Lamda[m][0]=A[m][m];
printf（"\n1Getalltheeigenvalues!
Congratulations!
\n"）;
fprintf（fp,"\n最大迭代次数k：
%4d\n",k）;
return（flag）;
}
elseif（m==0）
{
printf（"\n0Getalltheeigenvalues!
Congratulations!
\n"）;
fprintf（fp,"\n最大迭代次数k：
%4d\n",k）;
return（flag）;
}
else
gotoLoop1;
}
s1[0]=A[m-1][m-1]+A[m][m];s1[1]=0;
s2[0]=A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1];s2[1]=0;
root（s1,s2,Lamda,m）;
if（m==2）
{
for（i=0;i<2;i++）
{
Lamda[m][i]=s1[i];
Lamda[m-1][i]=s2[i];
}
printf（"\n2Getalltheeigenvalues!
Congratulations!
\n"）;
fprintf（fp,"\n最大迭代次数k：
%4d\n",k）;
return（flag）;
}
if（fabs（A[m-1][m-2]）<=Epsilon）
{
for（i=0;i<2;i++）
{
Lamda[m][i]=s1[i];
Lamda[m-1][i]=s2[i];
}
m=m-2;
gotoLoop2;
}
if（k==L）
{
printf（"\nLFailtogetalltheeigenvalues!
Pleasecheck!
\n"）;
flag=0;
return（flag）;
}
/*计算中间过度矩阵Mk*/
s=A[m-1][m-1]+A[m][m];
t=A[m-1][m-1]*A[m][m]-A[m-1][m]*A
升级会员 