数值分析大作业幂法反幂法特征值

《数值分析B 》大作业一

一． 算法设计方案：

先把带状线性矩阵A[501][501]转存为一个矩阵c[5][501].

1.1 λ1和λ501计算

1、判断λλ1501≠

首先使用幂法迭代求出矩阵的按模最大的特征值

0λ：若00>λ，则其必为max λ；于是可知max 501λλ=，接下来使用幂法迭代求出矩阵I A m ax λ-的按模最大的特征值：由矩阵理论可知I A max λ-的特征值按从小到大排列应为max 1λλ-、m ax 2λλ-、……、0。因此I A m ax λ-的按模最大的特征值应为m ax 1λλ-。又因为501max λλ=的值已求得，由此可求出1λ。

1.2

s λ的计算 s λ为矩阵A 按模最小的特征值，可以通过反幂法求得。

1.3 ik

λ的计算 ik λ可以对矩阵A 进行ik λ平移后，再用反幂法求出按模最小的特征值

λmin ，ik λ=

ik λ+λmin 。 1()||2n cond A λ

λ=，其中1λ和n λ分别是矩阵A 的模最大和最小特征值。 1.4det()A 的计算

矩阵A 的行列式可先对矩阵A 进行LU 分解后，det()A 等于U 所有对角线上元素的乘积。

二． 源程序

#include

#include

#include

#include

#include

#include

/*****************全局变量、函数声明*************/

#define N 501

#define EPSI 1.0e-12

#define r 2

#define s 2

double c[5][N]; /*A非零元素的压缩存储矩阵*/

double fuzhi(); /*对A赋值函数*/

void LUDet(); /*利用LU分解求解矩阵A的行列式*/

int int_max2(int a,int b); /*求两个数字中最大值的*/

int int_min2(int a,int b); /*求两个数字中最小值的*/

int int_max3(int a,int b,int c);/*求三个数字中最大值的*/

double mifa(); /*幂法计算矩阵A按模最大的特征值*/ double fmifa(); /*反幂法求矩阵A按模最大的特征值*/

//*主程序*//

void main()

{int i;

/*利用幂法计算矩阵A的最大特征值和最小特征值*/

double a1,a2;

fuzhi();

a1=mifa();

if (a1<0)

cout<<"矩阵A最小的特征值lambda1："<

else if(a1>=0)

cout<<"矩阵A最大的特征值lambda501："<

cout<

setprecision(12)<

for(i=0;i

for(i=0;i

c[2][i]=c[2][i]-a1;

a2=mifa()+a1;

if (a2<0)

cout<<"矩阵A最小的特征值lambda1："<

else if(a2>=0)

cout<<"矩阵A最大的特征值lambda501："<

cout<

setprecision(12)<

/*利用反幂法计算矩阵A的按模最小特征值*/

double a3;

fuzhi();

a3=fmifa();

cout<<"矩阵A按模最小的特征值lambdas:"<

cout<

setprecision(12)<

/*计算最接近特征值*/

fuzhi();

double u[39]={0};

cout<<"与数uk最接近的特征值："<

for(i=0;i<39;i++)

{

u[i]=a1+(i+1)*(a2-a1)/40;

c[2][i]=c[2][i]-u[i];

u[i]=fmifa()+u[i];

cout<<"lambda"<<

"[ik]"<<"("<

cout<

}

/*计算矩阵A的条件数，取2范数*/

double cond_A;

if(a1>0)

cond_A=fabs(a1/a3);

else if(a1<0)

cond_A=fabs(a2/a3);

cout<<"矩阵A的条件数:"<

cout<

setprecision(12)<

/*利用LU分解计算矩阵A的行列式*/

double det=1.0;

fuzhi();

LUDet(); //对矩阵A作LU分解

for(i=0;i

det*=c[2][i];

cout<<"矩阵A的行列式值："<

cout<

}

int int_max2(int a,int b)

{

return(a>b?a:b); /*求两个数字中最大值的子程序*/ }

int int_min2(int a,int b)

{

return(a

}

int int_max3(int a,int b,int c)

{ int t;

if(a>b)

t=a; /*求三个数字中最大值的子程序*/ else t=b;

if(t

return(t);

}

/*幂法求最大特征值*/

double mifa()

{

double u[N],unitu[N];

double sum,lengthu,beta,temp;

int i,j;

beta=0.0;

for(i=0;i

u[i]=0.5; /*迭代初始向量*/

do

{

temp=beta;

sum=0;

for(i=0;i

sum=sum+u[i]*u[i];

lengthu=sqrt(sum);

for(i=0;i

unitu[i]=u[i]/lengthu;

for(i=1;i<=N;i++)

{

sum=0;

for(j=int_max2(i-2,1);j<=int_min2(i+2,N);j++)

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

u[i-1]=sum;

}

sum=0;

for(i=0;i

sum=sum+unitu[i]*u[i];

beta=sum;

}while(fabs(beta-temp)>=EPSI);

return(beta);

}

/*反幂法求最小特征值*/

double fmifa()

{

double u[N],unitu[N],unitu_u[N];

double sum,lengthu,beta,temp;

int i,j;

beta=0.0;

for(i=0;i

u[i]=0.5; /*迭代初始向量*/ do

{

fuzhi();

temp=beta;

sum=0;

for(i=0;i

sum=sum+u[i]*u[i];

lengthu=sqrt(sum);

for(i=0;i

{

unitu[i]=u[i]/lengthu;

unitu_u[i]=unitu[i];

}

LUDet();

for(i=2;i<=N;i++)

{

sum=0;

for(j=int_max2(i-2,1);j<=i-1;j++)

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

unitu_u[i-1]-=sum;

}

u[N-1]=unitu_u[N-1]/c[2][N-1];

for(i=N-1;i>=1;i--)

{

sum=0;

for(j=i+1;j<=int_min2(i+2,N);j++)

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

u[i-1]=(unitu_u[i-1]-sum)/c[2][i-1];

}

sum=0;

for(i=0;i

sum=sum+unitu[i]*u[i];

beta=sum;

}while(fabs(beta-temp)>=EPSI);

beta=1/beta;

return(beta);

}

double fuzhi()

{

int i,j;

for(i=0;i<5;i++)

for(j=0;j

{

if(i==0||i==4)

{

c[i][j]=-0.064;

}

else if (i==1||i==3)

{

c[i][j]=0.16;

}

else if(i==2)

{

c[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))

-0.64*exp(0.1/(j+1));

}

}

return(c[i][j]);

}

/*LU三角分解*/

void LUDet()

{

double sum;

int k,i,j;

for(k=1;k<=N;k++)

{ for(j=k;j<=int_min2(k+s,N);j++)

{

sum=0;

for(i=int_max3(1,k-r,j-s);i<=k-1;i++)

sum+=c[k-i+s][i-1]*c[i-j+s][j-1];

c[k-j+s][j-1]-=sum;

}

for(j=k+1;j<=int_min2(k+r,N);j++)

{

sum=0;

for(i=int_max3(1,j-r,k-s);i<=k-1;i++)

sum+=c[j-i+s][i-1]*c[i-k+s][k-1];

c[j-k+s][k-1]=(c[j-k+s][k-1]-sum)/c[s][k-1];

}

}

}

三．程序结果：

四．分析初始向量选择对计算结果的影响

矩阵的初始向量选择，对结果的影响很大，选择不同的初始向量可能会得到不同阶的特征值。以幂法为例（反幂法原理相同），常见的初始向量选择有两种：

1、u[i]=1;(i=1,2,3...n)

2、u=0,u[i]=1(i为1-n任意整数)

若选择第1种初始向量，得到的λ1为：

若选择第2种初始向量，得到的λ1为：（i=0时）