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

数值分析之幂法及反幂法C语言程序实例1、算法设计方案：①求1、501和s的值：s：s表示矩阵的按模最小特征值，为求得s直接对待求矩阵A应用反幂法即可。1、501：已知矩阵A的特征值满足关系1n，要求1、及501时，可按如下方法求解：a．对矩阵A用幂法，求得按模最大的特征值1m。b．按平移量1m对矩阵A进行原点平移得矩阵1mBAI，对矩阵B用反幂法求得B的按模最小特征值2m。c．321mmm则：113min(,)mm，13max(,)nmm即为所求。②求和A的与数5011140kk最接近的特征值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和行列式detA：在（1）中用反幂法求矩阵A的按模最小特征值时，要用到Doolittle分解方法，在Doolittle分解完成后得到的两个矩阵分别为L和U，则A的行列式可由U阵求出，即：det(A)=det(U)。求得det(A)不为0，因此A为非奇异的实对称矩阵，则：max2()scondA，max和s分别为模最大特征值与模最小特征值。2、程序源代码：#includestdio.h#includestdio.h#includemath.h#defineN501//列#defineM5//行#defineR2//下带宽#defineS2//上带宽#defineK39#definee1.0e-12//误差限floatA[M][N];//初始矩阵floatu[N];//初始向量floaty[N],yy[N];floatmaximum,value1,value2,value_1,value_N,value_s,value_abs_max;constfloatb=0.16f,c=-0.064f;intmax_sign,max_position;voidInit_matrix_A()//初始化矩阵A{inti;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;}}voidInit_u()//初始化迭代向量{inti;for(i=0;iN;i++)u[i]=1.0;}voidGet_max()//获得绝对值最大的数值的模{inti;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;elsemax_sign=1;}voidGet_y()//单位化迭代向量{inti;for(i=0;iN;i++)y[i]=u[i]/maximum;}voidGet_u()//获得新迭代向量{inti;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];}voidGet_value()//获得迭代后特征值{value2=value1;value1=max_sign*u[max_position];}voidCheck_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;}}voidThe_value()//获取绝对值最大的特征值λ_501{Check_value();value_abs_max=value1;}voidThe_Other_value()//获取特征值λ_1{inti;floatvalue_temp=value1;for(i=0;iN;i++){A[2][i]-=value_temp;}Check_value();value1+=value_temp;if(value1value_temp){value_1=value1;value_N=value_temp;}else{value_N=value1;value_1=value_temp;}}intmin(inta,intb)//两值中取最小{if(ab)returna;elsereturnb;}intmax(inta,intb)//两值中取最大{if(ab)returnb;elsereturna;}voidResolve_LU(){intk,i,j,t;floattemp;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];}}}voidBack_substitution()//方程组回代过程{inti,t;floattemp=0;for(i=2;iN+1;i++){for(t=max(1,i-R);ti;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;i0;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];}}doubleDet_matrix()//求矩阵行列式值{inti;doubledet=1;Init_matrix_A();Resolve_LU();for(i=0;iN;i++)det=det*A[2][i];returndet;}floatGet_norm()//获得迭代向量模{inti;floatnormal=0;for(i=0;iN;i++)normal+=u[i]*u[i];normal=sqrt(normal);returnnormal;}voidGet_yy(floatnormal)//迭代向量单位化{inti;for(i=0;iN;i++){y[i]=u[i]/normal;yy[i]=y[i];}}voidGet_value_s()//获得绝对值最小的特征值{inti;value2=value1;value1=0;for(i=0;iN;i++)value1+=yy[i]*u[i];value1=1/value1;}voidValue_min()//反幂法求绝对值最小的特征值{floatnorm=0;intcount=0;value1=0,value2=0;Init_u();norm=Get_norm();Get_yy(norm);Back_substitution();Get_value_s();while(count10000){count++;norm=Get_norm();Get_yy(norm);Back_substitution();Get_value_s();if(fabs((value2-value1)/value1)e)break;}value_s=value1;}floatGet_cond_A()//求矩阵条件数{floatcond1;cond1=fabs(value_abs_max/value_s);returncond1;}voidValue_translation_min()//偏移条件下反幂法求特征值{inti,k;floattr;for(k=1;kK+1;k++){tr=value_1+k*(value_N-value_1)/40;Init_matrix_A();for(i=0;iN;i++)A[2][i]-=tr;Resolve_LU();Value_min();value_s+=tr;printf(k=%d=λi%d=%.13e\n,k,k,value_s);}}voidmain(){floatcond;doublevalue_det;printf(Contactme:[e-mail]\n);Init_matrix_A();//初始化矩阵AThe_value();//获取绝对值最大的特征值λ_501The_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的具有某些特征的特征值，主要用到的方法是幂法和反幂法，这两种方法从原理上看都是迭代法，因此迭代初始向量的选择对计算结果会产生一定影响，主要表现在收敛速度上。通过实际调试发现，对某些特殊的迭代初始值，确实对收敛结果及收敛速度产生影响，具体如下所列：以下结论建立在float数据类型基础之上；1.迭代初始值u[i]=c(i=1,2,…,501)且c的绝对值值极大（例如1.0e12以上），收敛结果可以稳定但收敛速度减慢，其原因为c的数量级与矩阵A中元素数量级差距过大，导致迭代次数以及运算量增大；2.迭代初始值u[i]=c(i=1,2,…,501)且c的绝对值值极小（例如1.0e-12以下），收敛结果并不稳定，且收敛速度减慢，其原因是计算机舍入误差将会影响计算结果；3.迭代初始值u[i](i=1,2,…,501)之间数量级偏差很大（例如1.0e12倍以上），收敛结果亦不稳定，且收敛速度减慢，其原因是人为使迭代过程中的权重发生较大区别，使迭代