北航数学学院数值分析作业1.docVIP

数值分析作业1 一、算法设计方案 计算λ1和λ501时，先用幂法求出A矩阵模最大的特征值，记为max1。然后对矩阵A 平移，求矩阵A-maxI 的模最大的特征值max2。若max20，则λ1=max1，λ501=max2+max1；否则，λ1=max1+max2；λ501=max1。 求λs时，用反幂法。由于反幂法需要解方程，故先对A三角分解（本例采取的是不列主元的Doolittle 分解），然后每解一次方程，就回代一次。 对于λik （k=1，2，···，39），利用带原点平移的反幂法，求矩阵A-μkI的模最小的特征值即可。 由于已经对A进行了三角分解，则求A的行列式det(A)时，只需求A所有对角元的乘积。其次因为A是非奇异的实对称矩阵，故A条件数Cond（A）2=│max1 / λs│。 二、源程序（Visual Studio 环境下的C语言） #includestdio.h #includestdlib.h #includemath.h #includemalloc.h #define N 501 #define s 2 /*对矩阵大小、 //#define s 2 带宽进行预定义，默认上下带宽一致 */ double Tmax(double (*p)[N+1]) ; /* 求矩阵的模最大的特征值，采用幂法 */ void Fjie( double (*p)[N+1] ); /** 对矩阵进行三角分解，采用不列主元 Doolittle分解法 **/ double *Hdai( double (*p)[N+1],double *q); /* 对已三角化的方程的系数矩阵求解，即为回代过程 */ double Tmin(double (*p)[N+1]) ; /* 求矩阵的模最小的特征值 ，采用反幂法 */ void main() { double a[6][502],b[6][502],x[40]; double t1,t2,max,min,det_a,cond_a,y; /************************************ t1--矩阵最小特征值； t2--矩阵最大特征值； max--模最大特征值； min--模最小特征值； det_a--a的行列式； cond_a--a的条件数； a存储的是题设矩阵A的非零元，a[i-j+s+1][j]=A[i][j]; ****************************************/ int i,j; for(i=1;i=502;i++) /* 对矩阵进行填充 */ { a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); if(i=3) {a[1][i]=-0.064;b[1][i]=-0.064;} if(i=2) { a[2][i]=0.16;b[2][i]=0.16;} if(i=500) { a[4][i]=0.16;b[4][i]=0.16;} if(i=499) {a[5][i]=-0.064;b[5][i]=-0.064;} } max=Tmax(a); /* 求A模最大的特征值 */ for(i=1;i=501;i++) b[3][i]=a[3][i]-max; t1=Tmax(b); if(max0) { t1+=max; t2=max; } else { t2=t1+max; t1=max; } /* 此时求得的t1，t2就是λ1，λ501 */ Fjie(a); for(i=1;i=39;i++) { y=t1+i*(t2-t1)/40; for(j=1;j=501;j++) b[3][j]=a[3][j]-y; x[i]=Tmin(b)+y; } /* x[i]的值就是 λii的值 */ for(i=1,det_a=1;i=501;i++) det_a=det_a*a[3][i]; min=Tmin(a); /* min的值就是λs的值 */ cond_a=fabs(max/min);