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

数值分析之幂法及反幂法C语言程序实例 1、算法设计方案： ①求、和的值： ：表示矩阵的按模最小特征值，为求得直接对待求矩阵A应用反幂法即可。 、：已知矩阵A的特征值满足关系 ，要求、及时，可按如下方法求解： 对矩阵A用幂法，求得按模最大的特征值。 按平移量对矩阵A进行原点平移得矩阵，对矩阵B用反幂法求得B的按模最小特征值。 则：，即为所求。 ②求和A的与数最接近的特征值（k=0，1，…39）： 求矩阵A的特征值中与最接近的特征值的大小，采用原点平移的方法： 先求矩阵 B=A-I 对应的按模最小特征值，则+即为矩阵A与最接近的特征值。 重复以上过程39次即可求得（k=0，1，…39）的值。 ③求A的（谱范数）条件数和行列式： 在（1）中用反幂法求矩阵A的按模最小特征值时，要用到Doolittle分解方法，在Doolittle分解完成后得到的两个矩阵分别为L和U，则A的行列式可由U阵求出，即：det(A)=det(U)。 求得det(A)不为0，因此A为非奇异的实对称矩阵，则： ，和分别为模最大特征值与模最小特征值。 2、程序源代码： #includestdio.h #includestdio.h #includemath.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;iN;i++) { A[0][i]= c; } for(i=1;iN;i++) { A[1][i]= b; } for(i=0;iN;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;iN-1;i++) { A[3][i]= b; } for(i=0;iN-2;i++) { A[4][i]= c; } } void Init_u() //初始化迭代向量 { int i; for(i=0;iN;i++) u[i]=1.0; } void Get_max() //获得绝对值最大的数值的模 { int i; max_position=0; maximum=fabs(u[0]); for(i=1;iN;i++) { if(maximumfabs(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;iN;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;iN-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() //幂法第二迭代格迭代