数值分析之幂法及反幂法C语言程序实例

数值分析之幂法及反幂法C 语言程序实例

1、算法设计方案：

①求1λ、501λ和s λ的值：

s λ：s λ表示矩阵的按模最小特征值，为求得s λ直接对待求矩阵A 应用反幂法即可。 1λ、501λ：已知矩阵A 的特征值满足关系 1n λλ<<，要求1λ、及501λ时，可

按如下方法求解： a ． 对矩阵A 用幂法，求得按模最大的特征值1m λ。

b ． 按平移量1m λ对矩阵A 进行原点平移得矩阵1m B

A I λ=+，对矩阵

B 用反幂法求得B 的按模最小特征值2m λ。

c ． 321m m m λλλ=-

则：113min(,)m m λλλ=，13max(,)n m m λλλ=即为所求。

②求和A 的与数5011

140k k λλμλ-=+最接近的特征值ik λ（k=0，1，…39）：

求矩阵A 的特征值中与k μ最接近的特征值的大小，采用原点平移的方法：

先求矩阵 B=A-k μI 对应的按模最小特征值k β，则k β+k μ即为矩阵A 与k μ最接近的特征值。

重复以上过程39次即可求得

ik λ（k=0，1，…39）的值。

③求A 的（谱范数）条件数2cond()A 和行列式det A ：

在（1）中用反幂法求矩阵A 的按模最小特征值时，要用到Doolittle 分解方法，在Doolittle 分解完成后得到的两个矩阵分别为L 和U ，则A 的行列式可由U 阵求出，即：det(A)=det(U)。 求得det(A)不为0，因此A 为非奇异的实对称矩阵，则： max 2()s

cond A λλ=，max λ和s λ分别为模最大特征值与模最小特征值。

2、程序源代码：

#include

#include

#include

#define N 501 //列

#define M 5 //行

#define R 2 //下带宽

#define S 2 //上带宽

#define K 39

#define e 1.0e-12 //误差限

float A[M][N]; //初始矩阵

float u[N]; //初始向量

float y[N],yy[N];

float maximum,value1,value2,value_1,value_N,value_s,value_abs_max;

const float b=0.16f,c=-0.064f;

int max_sign,max_position;

void Init_matrix_A() //初始化矩阵A

{

int i;

for(i=2;i

{

A[0][i]= c;

}

for(i=1;i

{

A[1][i]= b;

}

for(i=0;i

{

A[2][i]= (1.64-0.024*(i+1))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1));

}

for(i=0;i

{

A[3][i]= b;

}

for(i=0;i

{

A[4][i]= c;

}

}

void Init_u() //初始化迭代向量

{

int i;

for(i=0;i

u[i]=1.0;

}

void Get_max() //获得绝对值最大的数值的模

{

int i;

max_position=0;

maximum=fabs(u[0]);

for(i=1;i

{

if(maximum

{

max_position=i;

maximum=fabs(u[i]);

}

}

if(u[max_position]<0)

max_sign=-1;

else max_sign=1;

}

void Get_y() //单位化迭代向量

{

int i;

for(i=0;i

y[i]=u[i]/maximum;

}

void Get_u() //获得新迭代向量

{

int i;

u[0]=A[2][0]*y[0]+A[1][1]*y[1]+A[0][2]*y[2];

u[1]=A[3][0]*y[0]+A[2][1]*y[1]+A[1][2]*y[2]+A[0][3]*y[3];

u[N-2]=A[4][N-4]*y[N-4]+A[3][N-3]*y[N-3]+A[2][N-2]*y[N-2]+A[1][N-1]*y[N-1];

u[N-1]=A[4][N-3]*y[N-3]+A[3][N-2]*y[N-2]+A[2][N-1]*y[N-1];

for(i=2;i

u[i]=A[4][i-2]*y[i-2]+A[3][i-1]*y[i-1]+A[2][i]*y[i]+A[1][i+1]*y[i+1]+A[0][i+2]*y[i+2]; }