北航数值分析计算实习题目一 幂法 反幂法 求矩阵特征值.doc

数值分析计算实习题目一 SY0905303 解立垚 1.算法设计方案 对本题，首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，即可确定最小特征值λ1和最大特征值λ501。 使用反幂法求λs，其中需要解线性方程组，此处采用LU分解法解带状方程组。 通过带有原点平移的反幂法求出与数μk最接近的特征值λik。 Cond(A)=abs(λmax/λmin) ，其中λmax和λmin是按模最大和最小特征值。 分解矩阵A得出的LU矩阵，其中U矩阵的主对角线元素之积即为detA。 由于A的元素零元素较多，为节省储存量，将A的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A中的元素。 2.全部源程序 #include stdio.h #include math.h void init_a();//初始化A double get_an_element(int,int);//取A中的元素函数 double powermethod(double);//原点平移的幂法 double inversepowermethod(double);//原点平移的反幂法 int presolve(double);//三角LU分解 int solve(double [],double []);//解方程组 int max(int,int); int min(int,int); double (*u)[502]=new double[502][502];//上三角U数组 double (*l)[502]=new double[502][502];//单位下三角L数组 double a[6][502];//矩阵A int main() { int i,k; double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det; double lambda[40]; init_a();//初始化A lambdat1=powermethod(0); lambdat2=powermethod(lambdat1); lambda1=lambdat1lambdat2?lambdat1:lambdat2; lambda501=lambdat1lambdat2?lambdat1:lambdat2; presolve(0); lambdas=inversepowermethod(0); det=1; for(i=1;i=501;i++) det=det*u[i][i]; for (k=1;k=39;k++) { mu[k]=lambda1+k*(lambda501-lambda1)/40; presolve(mu[k]); lambda[k]=inversepowermethod(mu[k]); } printf(-----------结果汇总-------------\n); printf(λ=%1.11e λ=%1.11e\n,lambda1,lambda501); printf(λs=%1.11e\n,lambdas); printf(cond(A)=%1.11e\n,fabs(lambdat1/lambdas)); printf(detA=%1.11e \n,det); for (k=1;k=39;k++) {printf(λi%d=%1.11e ,k,lambda[k]);if(k % 3==0) printf(\n);} delete []u; delete []l;//释放堆内存 return 0; } void init_a()//初始化A { int i; for (i=3;i=501;i++) a[1][i]=a[5][502-i]=-0.064; for (i=2;i=501;i++) a[2][i]=a[4][502-i]=0.16; for (i=1;i=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); } double get_an_element(int i,int j)//从A中节省存储量的提取元素方法 { if (abs(i-j)=2) return a[i-j+3][j];else return 0; } double powermethod(double offset)//幂法 { int i,x1; double u[502],y[502]; double beta=0,prebeta=-1000,yita=0;