北航数值分析大作业一.docVIP

《数值分析》大作业一 1、算法设计方案： ⑴矩阵A的存储与检索： 将矩阵A[501][501],转存为新矩阵A[5][501]。在存入时，按每一行依次存入，存入时原矩阵与新矩阵的列号保持不变，其中原矩阵的主对角线上的元素存为新矩阵的第三行，最后所有元素存入新矩阵中。 ⑵求λ1，λ501，λs λs：λs表示矩阵的按模最小特征值，为求得λs直接对待求矩阵A应用反幂法即可。 λ1、λ501：已知矩阵A的特征值满足关系λ1＜…＜…λ501，要求λ1、及λ501时，可按如下方法求解： 对矩阵A用幂法，求得按模最大的特征值λn1。 按平移量λn1对矩阵A进行原点平移得矩阵，对矩阵B=A+λn1I用反幂法求得B的按模最小特征值λn2。 λn3=λn2一λn1 则：λ1=min（λn1，λn3），λ501=max（λn1，λn3） ⑶求A中与数μk=λ1+k（λ501-λ1）/40最接近的特征值λik（k=1，…39）： 先求矩阵 B=A-I 对应的按模最小特征值，则+即为矩阵A与最接近的特征值。 重复以上过程39次即可求得（k=0，1，…39）的值。 ⑷求A的（谱范数）条件数和行列式： 在（1）中用反幂法求矩阵A的按模最小特征值时，要用到Doolittle分解方法，在Doolittle分解完成后得到的两个矩阵分别为L和U，detA等于U所有对角线上元素的乘积。 ，λmax和λmin分别为模最大特征值与模最小特征值。 2、程序源代码： #includestdio.h #includemath.h #includestdlib.h #define N 501 /*定义列数*/ #define M 5 /*定义行数*/ #define e 1.0e-12 /*误差限*/ double A[M][N]; /*初始矩阵*/ double u[N]; /*初始向量*/ double y[N],yy[N]; double maximum,TZ1,TZ2,TZ_1,TZ_N,TZ_s,TZ_abs_max; int max_sign,max_position; /*对矩阵A进行初始化程序，存为A[5][501]*/ void CSH_A( ) { double b=0.16,c=-0.064; 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 CSH_u( ) { int i; for(i=0;iN;i++) u[i]=double (15*rand( )/32767); } /*幂法中对迭代向量进行单位化程序*/ void Get_y( ) { int i; for(i=0;iN;i++) y[i]=u[i]/maximum; } /*幂法中求解迭代向量的无穷范数*/ 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_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[