北航研究生数值分析B作业.docxVIP

数值分析计算实习题目一 设有的矩阵其中 矩阵的特征值满足 ,试求:的值的与数最接近的特征值的(谱范数)条件数和行列式.算法的设计方案 本题的核心算法是幂法、带原点平移的幂法、反幂法和LU分解法，要点在于选择算法时，应使的所有零元素都不存储。故算法设计的思路如下，第一步，对使用幂法(Powermethod)，可得的按模最大的特征值，记为；第二步，对使用带有原点平移的幂法，令平移量，可得另一端点的特征值记为；第三步，比较与的大小，根据条件可知， ；第四步，对使用反幂法(Inversepowermethod)，可得的按模最小的特征值（使用LU杜立特尔分解法）第五步，根据计算出,然后利用带有原点平移的反幂法求得,其中平移量,反幂法运算39次,可得;第六步, 根据定义,非奇异的实对称矩阵的谱范数条件数,其中分别是矩阵的模为最大和模为最小的特征值，对于本题，则有；第七步，由LU分解可知，，可得。由于题目要求算法中所有零元素均不为0，故构造一个givevalue（）函数为从一个的数组里为赋值。全部的源程序#include stdio.h#include math.hvoid init_a();//初始化Adouble givevalue(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];//矩阵Aint 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);}}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 givevalue(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;for (i=1;i=501;i++)u[i]=1,y[i]=0;//给向量u赋初值for (int