北航数值分析第一次大作业幂法反幂法.docx

北航数值分析第一次大作业幂法反幂法.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