北航数值分析第一次大作业幂法反幂法.docx
- 文档编号:12278632
- 上传时间:2023-04-17
- 格式:DOCX
- 页数:29
- 大小:68.11KB
北航数值分析第一次大作业幂法反幂法.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;i y[i]=u[i]/fabs(max); for(i=0;i { u[i]=0; for(j=0;j u[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;i yt[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;i ita=ita+pow(u[i],2); ita=sqrt(ita); for(i=1;i y[i]=u[i]/ita; LU(a,y,u);//调用LU求解线性方程组函数 if(k>0)//从k=1开始可以计算特征值 { beta=0; for(i=1;i beta=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;j a[i][j]=c[i][j]; } for(i=1;i b[i]=d[i]; for(k=1;k { for(i=k;i { s[i]=a[i][k]; for(t=1;t s[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;t u[k][j]=u[k][j]-l[k][t]*u[t][j]; } for(i=k+1;k l[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;t x[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;i for(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;j printf("%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;i printf("%f+i*%f\n",re[i],im[i]); printf("\nQR分解后的矩阵: \n"); for(i=0;i { for(j=0;j printf("%f",a[i][j]); printf("\r\n"); } //求特征向量 r=0; for(i=0
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 第一次 作业 幂法反幂法
![提示](https://static.bdocx.com/images/bang_tan.gif)