北航数值分析第一次大作业幂法反幂法.docx
《北航数值分析第一次大作业幂法反幂法.docx》由会员分享,可在线阅读,更多相关《北航数值分析第一次大作业幂法反幂法.docx(29页珍藏版)》请在冰豆网上搜索。
北航数值分析第一次大作业幂法反幂法
一、问题分析及算法描述
1.问题的提出:
(1)用幂法、反幂法求矩阵
的按摸最大和最小特征值,并求出相应的特征向量。
其中
要求:
迭代精度达到
。
(2)用带双步位移的QR法求上述的全部特征值,并求出每一个实特征值相应的特征向量。
2.算法的描述:
(1)幂法
幂法主要用于计算矩阵的按摸为最大的特征值和相应的特征向量。
其迭代格式为:
终止迭代的控制选用
。
幂法的使用条件为
实矩阵A具有n个线性无关的特征向量
,其相应的特征值
满足不等式
或
幂法收敛速度与比值
或
有关,比值越小,收敛速度越快。
(2)反幂法
反幂法用于计算
实矩阵A按摸最小的特征值,其迭代格式为:
每迭代一次都要求解一次线性方程组
。
当k足够大时,
,
可近似的作为矩阵A的属于
的特征向量。
比值
越小,收敛的越快。
反幂法要求矩阵A非奇异。
(3)带双步位移的QR分解法
QR方法适用于计算一般实矩阵的全部特征值,尤其适用于计算中小型实矩阵的全部特征值。
本算例中采用带双步位移的QR方法,可加速收敛,其迭代格式为:
二、计算结果及分析
1.计算结果:
(1)幂法:
初始条件:
最大迭代次数L=1000;向量
计算结果:
第1次迭代结果:
最大特征值:
0.00000e+000
第2次迭代结果:
最大特征值:
2.48910e+000相对误差:
1.00000e+000
第3次迭代结果:
最大特征值:
1.67719e+000相对误差:
4.84085e-001
第4次迭代结果:
最大特征值:
-2.10960e+000相对误差:
1.79503e+000
第5次迭代结果:
最大特征值:
-6.13203e-001相对误差:
2.44030e+000
…
…
第794次迭代结果:
最大特征值:
-1.97638e+000相对误差:
4.30074e-011
第795次迭代结果:
最大特征值:
-1.97638e+000相对误差:
3.04354e-013
********************最终迭代结果***************
特征值:
-1.97638e+000相对误差:
3.04354e-013
迭代次数:
795
(2)反幂法:
初始条件:
最大迭代次数L=1000;向量
运行结果:
第1次迭代结果:
最大特征值:
1.07542e+000
第2次迭代结果:
最大特征值:
-3.66550e+000相对误差:
1.29339e+000
第3次迭代结果:
最大特征值:
1.22709e+001相对误差:
1.29871e+000
第4次迭代结果:
最大特征值:
-1.03421e+000相对误差:
1.28650e+001
第5次迭代结果:
最大特征值:
-5.46339e-001相对误差:
8.92983e-001
…
…
第995次迭代结果:
最大特征值:
-3.24922e-001相对误差:
6.61387e-003
第996次迭代结果:
最大特征值:
-3.27176e-001相对误差:
6.88964e-003
第997次迭代结果:
最大特征值:
-3.29570e-001相对误差:
7.26527e-003
第998次迭代结果:
最大特征值:
-3.32147e-001相对误差:
7.75893e-003
第999次迭代结果:
最大特征值:
-3.34960e-001相对误差:
8.39573e-003
第1000次迭代结果:
最大特征值:
-3.38073e-001相对误差:
9.20985e-003
******************************
超过最大设定迭代次数,迭代失败!
(3)带双步位移的QR法:
初始条件:
最大迭代次数L=1000;向量
运行结果:
全部特征值:
0.747738+i*0.000000
-0.316674+i*0.025741
-0.316674+i*-0.025741
0.488580+i*0.139376
0.488580+i*-0.139376
1.045898+i*0.000000
-0.630981+i*0.000000
0.071427+i*0.539827
0.071427+i*-0.539827
1.265389+i*0.000000
-1.459878+i*0.000000
-1.521321+i*0.000000
-1.412499+i*0.148349
-1.412499+i*-0.148349
-1.976376+i*0.000000
1.810854+i*0.000000
1.362562+i*1.066117
1.362562+i*-1.066117
0.169262+i*1.909424
0.169262+i*-1.909424
特征向量(经谱范数归一化):
实特征值0.747738对应特征向量:
-0.062705-0.0223680.3043720.0644660.521833-0.1570240.136942-0.2181080.250264-0.043064
-0.228688-0.184632-0.0728710.1247210.0290700.102566-0.1363580.1677270.0857470.546165
实特征值1.045898对应特征向量:
-0.0180010.0196520.2734470.0705280.274896-0.1440150.0483850.376439-0.583051-0.054008
-0.168682-0.113430-0.0347090.0092040.4722910.125664-0.1906170.1131450.0462780.059871
实特征值-0.630981对应特征向量:
0.1068610.087709-0.024967-0.0208970.0643020.0340470.5351430.0463830.0288320.003479
-0.097276-0.3838010.089445-0.039560-0.036928-0.0213300.0148110.705836-0.1089040.082022
实特征值1.265389对应特征向量:
-0.0552010.0033990.2421910.1028470.372470-0.3728260.1139530.240659-0.310401-0.076590
-0.244632-0.192549-0.0772590.2633280.2016620.154166-0.4078140.1867820.0946490.173302
实特征值-1.459878对应特征向量:
0.427828-0.5468010.007822-0.3825800.0251990.0127880.0332410.005389-0.0040650.043524
-0.032112-0.0442330.135395-0.0065640.0012140.0201650.0116780.050001-0.5857650.013115
实特征值-1.521321对应特征向量:
0.236032-0.139250-0.0081430.638527-0.009049-0.002911-0.0013070.0030540.006515-0.030134
0.0127120.011368-0.018792-0.001753-0.005749-0.014290-0.005292-0.0145910.7175900.001369
实特征值-1.976376对应特征向量:
-0.227404-0.0481540.0226150.2973050.0703720.0399270.0785030.015822-0.0121820.605334
-0.083616-0.106270-0.573963-0.0199070.0038390.0513620.0365670.1156130.3327070.036954
实特征值1.810854对应特征向量:
-0.027768-0.051081-0.159642-0.054573-0.0844410.1183780.0295530.2110880.2038670.048627
0.0754700.026824-0.011103-0.584846-0.196009-0.0976680.673159-0.0368330.0033940.081244
2.结果分析
以上三种方法中,幂法计算共进行了795次迭代才达到收敛,计算量较大,收敛性不好;反幂法计算结果未能收敛,通过进一步分析发现,这是因为反幂法迭代程序未考虑按模最小特征值为复数的情况,造成迭代失败。
按双步位移的QR法迭代程序则考虑了特征值为复数的情况,并一次性给出了全部特征值。
三、结论
幂法、反幂法与QR法相比存在以下不足:
第一,需要设定初始条件,不适合自动运算;第二,幂法反幂法的迭代是否收敛依赖于特征值的分布情况,因此实际使用起来很不方便。
QR方法的收敛速度与特征值分布无关,是较好的特征值求解数值方法。
四、附录(VC++程序代码)
(1)幂法:
#include
#include
#defineN20
#definee1e-12
voidmain()
{
inti,j,k,r,rt;
doubleb,max,beta,betat,er,a[N][N],y[N],yt[N];
doubleu[N]={1};
printf("请输入最大迭代次数:
\n");
scanf("%d",&L);
for(i=0;i{
for(j=0;j{
if(i==j)a[i][j]=1.5*cos(double(i+1)+1.2*double(j+1));
elsea[i][j]=sin(0.5*double(i+1)+0.2*double(j+1));
}
}
for(k=0;;k++)
{
max=u[0];r=1;
for(i=1;i{
if(fabs(max){
max=u[i];r=i;
}
}
for(i=0;iy[i]=u[i]/fabs(max);
for(i=0;i{
u[i]=0;
for(j=0;ju[i]=u[i]+a[i][j]*y[j];
}
if(k>0)//从k=1开始可以计算特征值
{
printf("第%d次迭代结果:
\n",k);
beta=max*yt[rt];
printf("最大特征值:
%1.5e\n",beta);
}
if(k>1)//从k=2开始可以计算相对误差
{
er=fabs((beta-betat)/beta);
printf("相对误差:
%1.5e\n",er);
if(er{
printf("********************最终迭代结果*********************\n");
printf("特征值:
%1.5e\n",beta);
printf("相对误差:
%1.5e\n",er);
printf("迭代次数:
%d\n",k);
break;
}
if(k==L)//判断达到最大迭代次数是特征值是否收敛
if(er>e)
{
printf("超过最大设定迭代次数,迭代失败!
\n");break;
}
}
for(i=0;iyt[i]=y[i];
rt=r;betat=beta;
}
}
(2)反幂法:
#include
#include
#defineN20
#definee1e-12
voidLU(doublec[][N+1],doubleb[],doublex[]);//LU求解线性方程组函数的声明
voidmain()
{
inti,j,k,L;
doubleer,ita,beta,l,temp,a[N+1][N+1],y[N+1];
doubleu[N+1]={0,1};
for(i=1;i{
for(j=1;j{
if(i==j)a[i][j]=1.5*cos(double(i)+1.2*double(j));
elsea[i][j]=sin(0.5*double(i)+0.2*double(j));
}
}
printf("请输入最大迭代次数:
\n");
scanf("%d",&L);
for(k=0;;k++)
{
ita=0;
for(i=1;iita=ita+pow(u[i],2);
ita=sqrt(ita);
for(i=1;iy[i]=u[i]/ita;
LU(a,y,u);//调用LU求解线性方程组函数
if(k>0)//从k=1开始可以计算特征值
{
beta=0;
for(i=1;ibeta=beta+y[i]*u[i];
l=1/beta;
printf("第%d次迭代结果:
\n",k);
printf("最大特征值:
%1.5e\r\n",l);
}
if(k>1)//从k=2开始可以计算相对误差
{
er=fabs((1/beta-1/temp)*beta);
printf("相对误差:
%1.5e\n",er);
if(er{
printf("*********************************最终迭代结果***********************************\n");
printf("特征值:
%1.5e\n",beta);
printf("相对误差:
%1.5e\n",er);
printf("迭代次数:
%d\n",k);
break;
}
if(k==L)//判断达到最大迭代次数是特征值是否收敛
if(er>e)
{
printf("超过最大设定迭代次数,迭代失败!
\n");break;
}
}
temp=beta;
}
}
//LU求解线性方程组函数的定义
voidLU(doublec[][N+1],doubled[],doublex[])
{
doublea[N+1][N+1],b[N+1],y[N+1],l[N+1][N+1],u[N+1][N+1],s[N+1],temp,max;
inti,j,k,t=0,M[N+1];
for(i=1;i{
for(j=1;ja[i][j]=c[i][j];
}
for(i=1;ib[i]=d[i];
for(k=1;k{
for(i=k;i{
s[i]=a[i][k];
for(t=1;ts[i]=s[i]-l[i][t]*u[t][k];
}
max=0;
for(i=k;i{
if(max{
max=fabs(s[i]);
M[k]=i;
}
}
if(M[k]!
=k)
{
for(t=1;t{temp=l[k][t];l[k][t]=l[M[k]][t];l[M[k]][t]=temp;}
for(t=k;t{temp=a[k][t];a[k][t]=a[M[k]][t];a[M[k]][t]=temp;}
temp=s[k];s[k]=s[M[k]];s[M[k]]=temp;
}
u[k][k]=s[k];
for(j=k+1;k{
u[k][j]=a[k][j];
for(t=1;tu[k][j]=u[k][j]-l[k][t]*u[t][j];
}
for(i=k+1;kl[i][k]=s[i]/u[k][k];
}
for(k=1;k{
t=M[k];
temp=b[k];b[k]=b[t];b[t]=temp;
}
y[1]=b[1];
for(i=2;i{
y[i]=b[i];
for(t=1;t
y[i]=y[i]-l[i][t]*y[t];
}
x[N]=y[N]/u[N][N];
for(i=N-1;i>0;i--)
{
x[i]=y[i];
for(t=i+1;tx[i]=x[i]-u[i][t]*x[t];
x[i]=x[i]/u[i][i];
}
}
(3)带双步位移的QR法:
#include
#include
#include
#defineN20//稀疏矩阵大小
#defineL100//最大迭代次数
#definee1e-12//精度水平
voidHessenberg(doublea[][N]);//拟上三角函数的声明
voidQR(doubleb[][N],doublea[][N],intm);//QR分解函数的声明
voidLU(doublea[][N-1],doubleb[],doublex[]);//Doolittle分解函数的声明
//主函数
voidmain(void)
{
inti,j,k,m,r,n;
doubledet,tri,pd,sum;
doublevec[N],d[N-1],re[N],im[N],lam[N];
doublea[N][N],b[N][N],bp[N][N],c[N-1][N-1];
//系数矩阵
for(i=0;ifor(j=0;j{
if(i==j)a[i][j]=bp[i][j]=1.5*cos(double(i+1)+1.2*double(j+1));
elsea[i][j]=bp[i][j]=sin(0.5*double(i+1)+0.2*double(j+1));
}
Hessenberg(a);//调用拟上三角函数
printf("拟上三角化后的矩阵:
\n");
for(i=0;i{
for(j=0;jprintf("%f",a[i][j]);
printf("\n");
}
k=1;m=N-1;r=0;
third:
if(fabs(a[m][m-1])<=e)
{
re[r]=a[m][m];
im[r]=0;
r++;
m--;
}
else
gotofifth;
fourth:
if(m==0)
{
re[r]=a[0][0];
im[r]=0;
r++;
gotoeleventh;
}
if(m<0)
gotoeleventh;
if(m>0)
gotothird;
fifth:
det=a[m-1][m-1]*a[m][m]-a[m-1][m]*a[m][m-1];
tri=a[m-1][m-1]+a[m][m];
if(m==1||fabs(a[m-1][m-2])<=e)
{
pd=pow(tri,2)-4*det;
if(pd<0)
{
re[r]=re[r+1]=tri/2;
im[r]=sqrt(-pd)/2;
im[r+1]=-im[r];
r+=2;
}
else
{
re[r]=(tri+sqrt(pd))/2;
re[r+1]=(tri-sqrt(pd))/2;
im[r+1]=im[r]=0;
r+=2;
}
if(m==1)
gotoeleventh;
if(fabs(a[m-1][m-2])<=e)
{
m-=2;
gotofourth;
}
}
if(k==L)
{
printf("超过给定最大迭代次数,迭代失败!
\n");
gotoend;
}
else
gotonineth;
nineth:
for(i=0;i<=m;i++)
{
for(j=0;j<=m;j++)
{
sum=0;
for(n=0;n<=m;n++)
sum+=a[i][n]*a[n][j];
b[i][j]=sum-tri*a[i][j];
}
}
for(i=0;i<=m;i++)
b[i][i]+=det;
QR(b,a,m);//调用QR分解函数
k++;
gotothird;
eleventh:
printf("\n全部特征值:
\n");
for(i=0;iprintf("%f+i*%f\n",re[i],im[i]);
printf("\nQR分解后的矩阵:
\n");
for(i=0;i{
for(j=0;jprintf("%f",a[i][j]);
printf("\r\n");
}
//求特征向量
r=0;
for(i=0