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

数值分析计算实习题目一设有501501⨯的矩阵123499500501..........a b c b a b c c b a b c A cb a bc cb a b cba ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦其中 ()()()0.11.640.024s i n 0.20.641,2,,501;0.16,0.0ii a i i e i b c =--=⋅⋅⋅==-矩阵A 的特征值()1,2,,501i i λ=⋅⋅⋅满足12501,λλλ<<⋅⋅⋅< 1501||min ||s i i λλ≤≤=, 试求:1.1501,,s λλλ的值2. A 的与数5011140k kλλμλ-=+最接近的特征值(1,2,,39)k i k λ=3. A 的(谱范数)条件数2()cond A 和行列式det A .1. 算法的设计方案本题的核心算法是幂法、带原点平移的幂法、反幂法和LU 分解法，要点在于选择算法时，应使A 的所有零元素都不存储。

故算法设计的思路如下，第一步，对A 使用幂法(Powermethod)，可得A 的按模最大的特征值，记为a λ；第二步，对A 使用带有原点平移的幂法，令平移量1a p λ=，可得另一端点的特征值记为b λ； 第三步，比较a λ与b λ的大小，根据条件12501,λλλ<<⋅⋅⋅<可知{}1m i n ,a b λλλ=，{}501max ,a b λλλ=；第四步，对A 使用反幂法(Inversepowermethod)，可得A 的按模最小的特征值s λ（使用LU 杜立特尔分解法）第五步，根据5011140k kλλμλ-=+计算出k μ,然后利用带有原点平移的反幂法求得k i λ,其中平移量k k p μ=,反幂法运算39次,可得2239,,,i i i λλλ;第六步, 根据定义,非奇异的实对称矩阵A 的谱范数条件数()12||ncond A λλ=,其中1n λλ和分别是矩阵A 的模为最大和模为最小的特征值，对于本题，则有{}15012|max ||,|||()||s cond A λλλ=； 第七步，由LU 分解可知，A LU =，可得5011det det()det()det()(1501)n iii A LU L U u i ====⋅=≤≤∏。

由于题目要求算法中所有零元素均不为0，故构造一个givevalue （）函数为A 从一个6502⨯的数组里为A 赋值。

2. 全部的源程序#include <stdio.h> #include <math.h>void 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];//矩阵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=lambdat1<lambdat2?lambdat1:lambdat2; lambda501=lambdat1>lambdat2?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 k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;for (x1=1;x1<=501;x1++){u[x1]=0;for (int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(givevalue(x1,x2)-offset):givevalue(x1,x2))*y[x2];}prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出偏移量,误差,迭代次数}return (beta+offset);}double inversepowermethod(double offset)//反幂法{int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0; //给u[]赋初值for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];beta=1/beta;if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出偏移量,误差,迭代次数}return (beta+offset);}int presolve(double offset)//三角LU分解{int i,k,j,t;double sum;for (k=1;k<=501;k++)for (j=1;j<=501;j++){u[k][j]=l[k][j]=0;if (k==j) l[k][j]=1;} //初始化LU矩阵for (k=1;k<=501;k++){for (j=k;j<=min(k+2,501);j++){sum=0;for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++)sum=sum+l[k][t]*u[t][j];u[k][j]=((k==j)?(givevalue(k,j)-offset):givevalue(k,j))-sum;}if (k==501) continue;for (i=k+1;i<=min(k+2,501);i++){sum=0;for (t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(givevalue(i,k)-offset):givevalue(i,k))-sum)/u[k][k];}}return 0;}int solve(double x[],double b[])//求解方程组{int i,t;double y[502];double sum;y[1]=b[1];for (i=2;i<=501;i++){sum=0;for (t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for (i=500;i>=1;i--){sum=0;for (t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}int max(int x,int y){return (x>y?x:y);}int min(int x,int y){return (x<y?x:y);}3. 计算结果计算结果如图1所示图1偏移量offset ，误差err ，迭代次数k 也在结果中有所反应,部分结果如图2所示图24. 讨论迭代初始向量的选取对计算结果的影响,并说明原因.在main 程序中,为初始向量u 赋初值[]1(1501)u i i =≤≤得到上述结果. 下面改变初始向量u 的赋值,观察计算结果的差异[][]1(1),0(2501)u i i u j j ===≤≤时,结果如图3图3[][]=≤≤=≤≤时,结果如图41(1101),0(102501)u i i u j j图4[][]u i i u j j=≤≤=≤≤时,结果如图51(1201),0(202501)图5[][]=≤≤=≤≤时,结果如图6u i i u j j1(1301),0(302501)图6[][]=≤≤=≤≤时,结果如图7u i i u j j1(1401),0(402501)图7[][]u i i u j j=≤≤==时,结果如图80(1500),1(501)图8通过观察发现,矩阵迭代法的初始向量选取非常重要,并不是热选的初始向量都能收敛得到想要的结果,不同的初始向量可能是计算结果收敛于不同接的特征值与特征向量,而不一定收敛于第一阶.{}4,3,2,1A diag =,其特征值显然为12344,3,2,1λλλλ====,相应的特征向量为1234(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),T T T T X X X X ====通过不同的初始向量得到的迭代表1从表1中可以看到:由于选取初始向量的不同,矩阵迭代的结果不一定收敛于X1与λ1,如表中②、⑤栏所示;若初始向量很接近于某一阶特征向量,如③、④、⑤、⑥栏中的情况,则收敛结果可能是最接近的这阶向量,也可能收敛于其他阶特征向量。