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

数值分析之幂法及反幂法C 语言程序实例1、算法设计方案：①求1λ、501λ和s λ的值：s λ：s λ表示矩阵的按模最小特征值，为求得s λ直接对待求矩阵A 应用反幂法即可。

1λ、501λ：已知矩阵A 的特征值满足关系 1n λλ<<，要求1λ、及501λ时，可按如下方法求解： a ． 对矩阵A 用幂法，求得按模最大的特征值1m λ。

b ． 按平移量1m λ对矩阵A 进行原点平移得矩阵1m BA I λ=+，对矩阵B 用反幂法求得B 的按模最小特征值2m λ。

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

②求和A 的与数5011140k 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()scond A λλ=，max λ和s λ分别为模最大特征值与模最小特征值。

2、程序源代码：#include<stdio.h>#include<stdio.h>#include<math.h>#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<N;i++){A[0][i]= c;}for(i=1;i<N;i++){A[1][i]= b;}for(i=0;i<N;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<N-1;i++){A[3][i]= b;}for(i=0;i<N-2;i++){A[4][i]= c;}}void Init_u() //初始化迭代向量{int i;for(i=0;i<N;i++)u[i]=1.0;}void Get_max() //获得绝对值最大的数值的模{int i;max_position=0;maximum=fabs(u[0]);for(i=1;i<N;i++){if(maximum<fabs(u[i])){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<N;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<N-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]; }void Get_value() //获得迭代后特征值{value2=value1;value1=max_sign*u[max_position];}void Check_value() //幂法第二迭代格迭代{Init_u();Get_max();Get_y();Get_u();Get_value();while(1){Get_max();Get_y();Get_u();Get_value();if(fabs((value2-value1)/value1)<e)break;}}void The_value() //获取绝对值最大的特征值λ_501 {Check_value();value_abs_max=value1;}void The_Other_value() //获取特征值λ_1{int i;float value_temp=value1;for(i=0;i<N;i++){A[2][i]-=value_temp;}Check_value();value1+=value_temp;if(value1<value_temp){value_1=value1;value_N=value_temp;}else{value_N=value1;value_1=value_temp;}}int min(int a,int b) //两值中取最小{if(a<b)return a;elsereturn b;}int max(int a,int b) //两值中取最大{if(a<b)return b;elsereturn a;}void Resolve_LU(){int k,i,j,t;float temp;for(k=1;k<=N;k++){for(j=k;j<=min(k+S,N);j++){temp=0;for(t=max(max(1,k-R),j-S);t<=k-1;t++)temp+=A[k-t+S][t-1]*A[t-j+S][j-1];A[k-j+S][j-1]=A[k-j+S][j-1]-temp;}for(i=k+1;i<=min(k+R,N);i++){temp=0;for(t=max(max(1,i-R),k-S);t<=k-1;t++)temp+=A[i-t+S][t-1]*A[t-k+S][k-1];A[i-k+S][k-1]=(A[i-k+S][k-1]-temp)/A[S][k-1];}}}void Back_substitution()//方程组回代过程{int i,t;float temp=0;for(i=2;i<N+1;i++){for(t=max(1,i-R);t<i;t++)y[i-1]-=A[i-t+S][t-1]*y[t-1];}u[N-1]=y[N-1]/A[S][N-1];for(i=N-1;i>0;i--){temp=0;for(t=i+1;t<=min(i+S,N);t++)temp+=A[i-t+S][t-1]*u[t-1];u[i-1]=(y[i-1]-temp)/A[S][i-1];}}double Det_matrix() //求矩阵行列式值{int i;double det=1;Init_matrix_A();Resolve_LU();for(i=0;i<N;i++)det=det*A[2][i];return det;}float Get_norm() //获得迭代向量模{int i;float normal=0;for(i=0;i<N;i++)normal+=u[i]*u[i];normal=sqrt(normal);return normal;}void Get_yy(float normal) //迭代向量单位化{int i;for(i=0;i<N;i++){y[i]=u[i]/normal;yy[i]=y[i];}}void Get_value_s() //获得绝对值最小的特征值{int i;value2=value1;value1=0;for(i=0;i<N;i++)value1+=yy[i]*u[i];value1=1/value1;}void Value_min() //反幂法求绝对值最小的特征值{float norm=0;int count=0;value1=0,value2=0;Init_u();norm=Get_norm();Get_yy(norm);Back_substitution();Get_value_s();while(count<10000){count++;norm=Get_norm();Get_yy(norm);Back_substitution();Get_value_s();if(fabs((value2-value1)/value1)<e)break;}value_s=value1;}float Get_cond_A() //求矩阵条件数{float cond1;cond1=fabs(value_abs_max/value_s);return cond1;}void Value_translation_min() //偏移条件下反幂法求特征值{int i,k;float tr;for(k=1;k<K+1;k++){tr=value_1+k*(value_N-value_1)/40;Init_matrix_A();for(i=0;i<N;i++)A[2][i]-=tr;Resolve_LU();Value_min();value_s+=tr;printf("k=%d =>>>λi%d=%.13e\n",k,k,value_s);}}void main(){float cond;double value_det;printf("Contactme:****************\n");Init_matrix_A(); //初始化矩阵AThe_value(); //获取绝对值最大的特征值λ_501 The_Other_value(); //获取特征值λ_1printf("λ1=%.13e\n",value_1);printf("λ501=%.13e\n",value_N);value_det=Det_matrix(); //求矩阵行列式值Value_min(); //反幂法求绝对值最小的特征值printf("λs=%.13e\n",value_s);cond=Get_cond_A(); //求矩阵条件数Value_translation_min();//偏移条件下反幂法求特征值printf("cond_A=%.13e\n",cond);printf("value_det=%.13e\n",value_det);}3、程序运行结果：4、迭代初始向量的选取对计算结果的影响：本次计算实习求矩阵A的具有某些特征的特征值，主要用到的方法是幂法和反幂法，这两种方法从原理上看都是迭代法，因此迭代初始向量的选择对计算结果会产生一定影响，主要表现在收敛速度上。