北航数值分析大作业2.docx
- 文档编号:8624000
- 上传时间:2023-02-01
- 格式:DOCX
- 页数:23
- 大小:24.59KB
北航数值分析大作业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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业