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

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

《北航数值分析大作业 第一题 幂法与反幂法.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业 第一题 幂法与反幂法.docx（18页珍藏版）》请在冰豆网上搜索。

北航数值分析大作业第一题幂法与反幂法

《数值分析》计算实习题目

-----第一题

学院：

学号：

姓名：

北京航空航天大学

2011年10月

一、算法设计方案

矩阵A是上半带宽为2、下半带宽为2的带状矩阵。

为节省内存空间，首先将矩阵

压缩到矩阵

中，在矩阵C中检索矩阵A的带内元素的方法是：

1.1求

和

的值

首先调用幂法函数求出矩阵A按模最大的特征值

，再以

为平移量再次调用幂法函数求出按模最大特征值

。

比较

和

的大小，较大值为

，较小值为

。

1.2求

的值

调用反幂法函数求出矩阵A按模最小特征值

。

1.3求与数

最接近的特征值

。

以

为平移量调用39次反幂法函数求出39个按模最小特征值

，则所求与

最接近的特征值

。

1.4求A的（谱范数）条件数

和行列式

。

由于A是非奇异的实对称矩阵，故

，其中

和

是矩阵A的按模最大和按模最小特征值，以由1.1和1.2求出；对矩阵A进行LU分解，A的行列式值即为U矩阵对角线元素的乘积。

二、算法源程序

#include

#include

#include

#include

doubleLamda_1,Lamda_501,Lamda_s,q,u,L[40],condA,detA;

inti,j,k,t,f;

classsolution//定义“解题”类

{

private:

//声明private型数据

intm[501];

doublea[5][501],c[5][501],v[501],b,sum,max,p,y[501],r,d,x[501];

public:

voidsolution_input（）;//声明数据输入函数

doublesolution_powermethod（doubler）;//声明幂法函数，r为平移量

doublesolution_inversepowermethod（doubler）;//声明反幂法函数，r为平移量

doublesolution_detA（）;//声明求行列式函数

doubles_min（doublea,doubleb）;//声明取最小值函数

doubles_max（doublea,doubleb）;//声明取最大值函数

~solution（）//在运算结束后调用析构函数释放内存

{

}

};

voidmain（）

{

solutionsolve;//创建solve对象

solve.solution_input（）;//调用数据输入函数给矩阵A赋值

Lamda_1=solve.solution_powermethod（0）;//以0为平移量，调用幂法函数

//以Lamda_1为平移量，调用幂法函数

Lamda_501=solve.solution_powermethod（Lamda_1）;

if（Lamda_1>Lamda_501）//两个特征值大的为Lamda_501,较小的为Lamda_1

{

q=Lamda_1;

Lamda_1=Lamda_501;

Lamda_501=q;

}

//调用反幂法函数，求得按模最小特征值

Lamda_s=solve.solution_inversepowermethod（0）;

cout<<"************第1小题结果*************"<

//结果分别输出到屏幕和txt文件中

cout<<"Lamda_1="<

:

scientific）<

cout<<"Lamda_501="<

:

scientific）<

cout<<"Lamda_s="<

:

scientific）<

ofstreamoutput（"Result.txt"）;

output<<"************第1小题结果*************"<

output<<"Lamda_1="<

:

scientific）<

output<<"Lamda_501="<

:

scientific）<

output<<"Lamda_s="<

:

scientific）<

cout<<"************第2小题结果*************"<

output<<"************第2小题结果*************"<

for（f=1;f<40;f++）//将u的值作为自变量，调用反幂法函数，求得与u最接近的特征值

{

u=Lamda_1+f*（Lamda_501-Lamda_1）/40;

L[f]=solve.solution_inversepowermethod（u）;

cout<<"L["<

:

scientific）<

output<<"L["<

:

scientific）<

}

detA=solve.solution_detA（）;//调用求行列式函数求得detA

cout<<"************第3小题结果*************"<

output<<"************第3小题结果*************"<

cout<<"detA="<

:

scientific）<

output<<"detA="<

:

scientific）<

if（detA!

=0）//若A为非奇异矩阵，求得A的条件数

{

condA=fabs（Lamda_1/Lamda_s）;

cout<<"condA="<

:

scientific）<

output<<"condA="<

:

scientific）<

}

cout<<"***************实习昨天第一题完成******************"<

output<<"***************实习昨天第一题完成******************"<

}

//定义输入数据函数

voidsolution:

:

solution_input（）

{

/**********将矩阵A存放到c[5][501]中*******/

for（i=0;i<5;i++）

{

for（j=0;j<501;j++）

{c[i][j]=0;}

}

for（j=0;j<501;j++）

{c[2][j]=（（1.64-0.024*（j+1））*sin（0.2*（j+1））-0.64*exp（0.1/（j+1）））;}

for（j=1;j<501;j++）

{c[1][j]=0.16;}

for（j=2;j<501;j++）

{c[0][j]=-0.064;}

for（j=0;j<500;j++）

{c[3][j]=0.16;}

for（j=0;j<499;j++）

{c[4][j]=-0.064;}

/**********将矩阵c[5][501]备份*******/

for（i=0;i<5;i++）

{

for（j=0;j<501;j++）

{a[i][j]=c[i][j];}

}

}

//定义幂法函数，自变量r是平移量

doublesolution:

:

solution_powermethod（doubler）

{

/**********初始化特征向量*******/

for（i=0;i<501;i++）

{v[i]=1;}

b=0;

/**********幂法迭代过程*******/

for（;;）

{

sum=0;

for（i=0;i<501;i++）

{sum+=pow（v[i],2）;}

p=sqrt（sum）;

for（i=0;i<501;i++）

{y[i]=v[i]/p;}

for（i=0;i<501;i++）

{

sum=0;

for（j=0;j<501;j++）

{

if（fabs（i-j）<3）

{

if（i==j）

sum+=（c[i-j+2][j]-r）*y[j];

else

sum+=c[i-j+2][j]*y[j];

}

}

v[i]=sum;

}

d=0;

for（i=0;i<501;i++）

{

d+=y[i]*v[i];

}

if（fabs（d-b）/fabs（d）<=1e-12）//设置迭代精度

break;

b=d;

}

returnd+r;//返回值即为按模最大特征值

}

//定义反幂法函数，自变量r为平移量

doublesolution:

:

solution_inversepowermethod（doubler）

{

/**********初始化各向量*******/

for（j=0;j<501;j++）

{x[j]=0;}

for（j=0;j<501;j++）

{y[j]=0;}

for（i=0;i<5;i++）

{

for（j=0;j<501;j++）

{c[i][j]=a[i][j];}

}

for（j=0;j<501;j++）

{

c[2][j]=c[2][j]-r;

}

for（i=0;i<501;i++）

{v[i]=1;}

b=0;

/**********先对矩阵A进行Doolittle分解，分解所得的L、U矩阵仍存放到c矩阵中*******/

for（k=0;k<501;k++）

{

for（j=k;j<=solution:

:

s_min（k+2,500）;j++）

{

sum=0;

for（t=solution:

:

s_max（0,solution:

:

s_max（k-2,j-2））;t<=k-1;t++）

{sum+=c[k-t+2][t]*c[t-j+2][j];}

c[k-j+2][j]=c[k-j+2][j]-sum;

}

for（i=k+1;i<=solution:

:

s_min（k+2,500）;i++）

{

if（k>=500）

break;

sum=0;

for（t=solution:

:

s_max（0,solution:

:

s_max（k-2,i-2））;t<=k-1;t++）

{sum+=c[i-t+2][t]*c[t-k+2][k];}

c[i-k+2][k]=（c[i-k+2][k]-sum）/c[2][k];

}

}

/**********求第k代特征向量y，并备份*******/

for（;;）

{

sum=0;

for（i=0;i<501;i++）

{sum+=pow（v[i],2）;}

q=sqrt（sum）;

for（i=0;i<501;i++）

{

y[i]=v[i]/q;

x[i]=y[i];

}

/**********用Doolittle分解法求解第k代迭代向量v以及特征值*******/

for（i=1;i<501;i++）

{

sum=0;

for（t=solution:

:

s_max（0,i-2）;t<=（i-1）;t++）

{sum+=c[i-t+2][t]*y[t];}

y[i]=y[i]-sum;

}

v[500]=y[500]/c[2][500];

for（i=499;i>=0;i--）

{

sum=0;

for（t=（i+1）;t<=solution:

:

s_min（i+2,500）;t++）

{sum+=c[i-t+2][t]*v[t];}

v[i]=（y[i]-sum）/c[2][i];

}

d=0;

for（i=0;i<501;i++）

{

d+=x[i]*v[i];

}

if（fabs（（1/d）-（1/b））/fabs（1/d）<=1e-12）//设置迭代精度

break;

b=d;

}

returnr+1/d;//返回值即为按模最小特征值

}

//定义求行列式函数

doublesolution:

:

solution_detA（）

{

/********将备份的矩阵A的值赋给矩阵c*********/

for（i=0;i<5;i++）

{

for（j=0;j<501;j++）

{c[i][j]=a[i][j];}

}

/**********先对矩阵A进行Doolittle分解，分解所得的L、U矩阵仍存放到c矩阵中*******/

for（k=0;k<501;k++）

{

for（j=k;j<=solution:

:

s_min（k+2,500）;j++）

{

sum=0;

for（t=solution:

:

s_max（0,solution:

:

s_max（k-2,j-2））;t<=k-1;t++）

{sum+=c[k-t+2][t]*c[t-j+2][j];}

c[k-j+2][j]=c[k-j+2][j]-sum;

}

for（i=k+1;i<=solution:

:

s_min（k+2,500）;i++）

{

if（k>=500）

break;

sum=0;

for（t=solution:

:

s_max（0,solution:

:

s_max（k-2,i-2））;t<=k-1;t++）

{sum+=c[i-t+2][t]*c[t-k+2][k];}

c[i-k+2][k]=（c[i-k+2][k]-sum）/c[2][k];

}

}

/**********detA为矩阵u对角线上运算的乘积*******/

sum=1;

for（j=0;j<501;j++）

{sum*=c[2][j];}

returnsum;

}

//定义取最小值函数

doublesolution:

:

s_min（doublea,doubleb）

{

if（a<=b）returna;

elsereturnb;

}

//定义取最大值函数

doublesolution:

:

s_max（doublea,doubleb）

{

if（a>=b）returna;

elsereturnb;

}

三、算法运行结果

将计算结果分别输出到屏幕和txt文件中，如下所示：

************第1小题结果*************

Lamda_1=-1.070011361502e+001

Lamda_501=9.724634098777e+000

Lamda_s=-5.557910794230e-003

************第2小题结果*************

L[1]=-1.018293403315e+001

L[2]=-9.585707425068e+000

L[3]=-9.172672423928e+000

L[4]=-8.652284007898e+000

L[5]=-8.093483808675e+000

L[6]=-7.659405407692e+000

L[7]=-7.119684648691e+000

L[8]=-6.611764339397e+000

L[9]=-6.066103226595e+000

L[10]=-5.585101052628e+000

L[11]=-5.114083529812e+000

L[12]=-4.578872176865e+000

L[13]=-4.096470926260e+000

L[14]=-3.554211215751e+000

L[15]=-3.041090018133e+000

L[16]=-2.533970311130e+000

L[17]=-2.003230769563e+000

L[18]=-1.503557611227e+000

L[19]=-9.935586060075e-001

L[20]=-4.870426738850e-001

L[21]=2.231736249575e-002

L[22]=5.324174742069e-001

L[23]=1.052898962693e+000

L[24]=1.589445881881e+000

L[25]=2.060330460274e+000

L[26]=2.558075597073e+000

L[27]=3.080240509307e+000

L[28]=3.613620867692e+000

L[29]=4.091378510451e+000

L[30]=4.603035378279e+000

L[31]=5.132924283898e+000

L[32]=5.594906348083e+000

L[33]=6.080933857027e+000

L[34]=6.680354092112e+000

L[35]=7.293877448127e+000

L[36]=7.717111714236e+000

L[37]=8.225220014050e+000

L[38]=8.648666065193e+000

L[39]=9.254200344575e+000

************第3小题结果*************

detA=2.772786141752e+118

condA=1.925204273902e+003

***************实习作业第一题完成******************

四、讨论迭代初始向量的选取对计算结果的影响

矩阵的初始向量选择，对结果的影响很大，选择不同的初始向量可能会得到不同阶的特征值。

以幂法为例（反幂法原理相同），常见的初始向量选择有两种：

1、

2、

，此处以

为例讨论。

通过第一种方法构造迭代初始向量，求得第一小题的结果是：

************第1小题结果*************

Lamda_1=-1.070011361502e+001

Lamda_501=9.724634098777e+000

Lamda_s=-5.557910794230e-003

通过第二种方法构造迭代初始向量，求得第二小题的结果是：

************第1小题结果*************

Lamda_1=-2.080981085336e+000

Lamda_501=9.978750038032e-001

Lamda_s=2.668886923785e-002

结果发现两种计算结果相差较大，而第一种初始向量得到的特征值模更大，实际上，此种初始化向量的方法能得到最准确的按模最大特征值，但迭代数次会较第2种方法多很多。

当增加迭代初始向量中非零元素的个数时，计算出的结果会越接近准确值。

实际上，当初始向量

中的0元素较多时，可能

的情况较为普遍，许多

都有可能等于0，此时计算出的结果便与最大特征值差距较大。

实际操作中可以选择不同形式的初始向量，利用求得的特征值，找到其中最大的即为需求的按模最大特征值。