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

数值分析之幕法及反幕法C语言程序实例1、算法设计方案:①求1、 501和s的值:S： s表示矩阵的按模最小特征值，为求得s直接对待求矩阵A应用反幕法即可。1、 501:已知矩阵A的特征值满足关系III，要求1、及501时，可按如下方法求解:a.对矩阵A用幕法，求得按模最大的特征值ml°b.按平移量 m1对矩阵A进行原点平移得矩阵m11，对矩阵B用反幕法c.则:求得B的按模最小特征值a.对矩阵A用幕法，求得按模最大的特征值ml°b.按平移量 m1对矩阵A进行原点平移得矩阵m11，对矩阵B用反幕法c.则:求得B的按模最小特征值m3m2mlm2°1min(m1,m3)nmax(m1,m3)即为所求。②求和A的与数计最接近的特征值ik(k=0.1，…39):求矩阵A的特征值中与k最接近的特征值的大小，采用原点平移的方法:先求矩阵B=A-kI对应的按模最小特征值k，则k+k即为矩阵A与k最接近的…39…39)的值。重复以上过程39次即可求得 ik(k=0,1,③求A的(谱数)条件数cond(A)2和行列式detA:在(1)中用反幕法求矩阵 A的按模最小特征值时，要用到 Doolittle分解方法，在Doolittle 分解完成后得到的两个矩阵分别为 L和U,贝UA的行列式可由U阵求出，即：det(A)=det(U)°求得det(A)不为0,因此A为非奇异的实对称矩阵，则：cond(A)2 吕,max和S分别为模最大特征值与模最小特征值。S2、程序源代码：#include<stdio.h>#include<stdio.h>#include<math.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;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;}}voidInit_u() //初始化迭代向量{inti;for(i=0;i<N;i++)u[i]=1.0;}voidGet_max() //获得绝对值最大的数值的模{inti;max_position=0;maximum=fabs(u[O]);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;elsemax_sign=1;}voidGet_y() //单位化迭代向量{inti;for(i=0;i<N;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;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];}voidGet_value() //获得迭代后特征值{value2=value1;value仁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;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;}}//两值中取最小////两值中取最小//两值中取最大{if(a<b)returna;elsereturnb;}intmax(inta,intb){if(a<b)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;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];}}doubleDet_matrix()// 求矩阵行列式值{inti;doubledet=1;Init_matrix_A();Resolve_LU();for(i=0;i<N;i++)det=det*A[2][i];returndet;}floatGet_norm() //获得迭代向量模{inti;floatnormal=0;for(i=0;i<N;i++)normal+=u[i]*u[i];normal=sqrt(normal);returnnormal;}voidGet_yy(floatnormal) //迭代向量单位化{inti;for(i=0;i<N;i++){y[i]=u[i]/normal;yy[i]=y[i];}}voidGet_value_s()//获得绝对值最小的特征值{inti;value2=value1;value1=0;for(i=0;i<N;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(count<10000){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;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);}}voidmain(){floatcond;doublevalue_det;printf("Contactme:731363227qq.\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、程序运行结果:ofNuaericalAnaly.口以*C:\DocuMentsandSettings\Ad>inistrator\^：面ITasik1ESH丄财1聊刘湄数值分析大徐if丄ofNuaericalAnaly.0730092315G?4e+0911501=9-7247S1445312Se+ao0ls=-5.55?9231±20646e-0&3k=l=»>Xil=-l.01829338B0482e+061=»> 12=-9,585707B?82456tr*0e0J<=3=>»ki3=-9.17267238B6095e+060k=4=>»>i4=-R.fiS22R3717R7R2pr+nPHk=5->»Xi5=-8.09?4837^e4035e+000kf=»>Xi6=-?*6594052538276e+060](=?=»>y.i7—=■>»X16—6h-9=»>Xi9=-&.0Gtl033093929c+000k-10-»>X .5851011B34101e+000k=ll=»>\=»>入il2=-4.5?88?22643629e+000k=13=»> il3=-4.B964710377157c+000k=14=»>X.il^=-3.5542112756521e+000kT5=»>入±15^-3・9410?aS27312?e*008k=1ft=>»Xil6=-2.52643ft3ft4721fie+HRWk-1?-»>Xil7^-2.003233fl47"il83e*000k=18w”〉入i!8=-l.5035575726070c+000](=19=>〉》入il9=-9-9355ebl73B874e-001k-20-»>Xi2e—4.8?042t746?838e-001h-21->>>X121-B_231736JS223G7c602k=22-»>X122-5-3241?47G17722e-001k=S3=»>'123=1-0S28989899904e+000k=24=»>入124=1.589445903S9?3e+000k=25=>»'、125=2.06033B4t35733e+008k=26=»>X126=2.5580755751580e+B00k=27=»>X .08O24a402B364e*S00=>»'、、i2R=3.k=2?-»>Xi29-4*09137S5751?0?e+P00k-30->»X130=4.60303552821280+600k=3t->»入131=5.13292419&3103e*600k-32-»>入132-5-594906423240?e也001^93->>>X133H；-0e0933?347?4Ge*600k=34-»>X_6803541369736e+000k=3S=»>h135=7.2938773855567e+000k=36=»>X136=7.7171115fcl4474e+800k=37=»>hi3V=8-225219?4i4637e+000k=38=>»人