北航数值分析大作业一(纯原创,高分版)

数值分析上机实习作业一题目：设有501×501的实对称矩阵A，A其中ai=1.64-0.024isin0.2λ求λ1,和λ501的值。求A的与数μk=λ求A的（普范数）条件数cond(A)和行列式detA。说明在所用的算法中，凡是要给出精度水平的ε，都取ε=10选择算法的时候应使矩阵A的所有零元素都不存储。打印以下内容：(1)算法设计方案和思路。(2)全部源程序。(3)特征值λ1,λ501,λs,λik(4)发现的现象与问题。算法设计思路和方案因为题目只要求求部分特征值，及最大的特征值和最小的特征值以及距离某些给定实数最近的特征值，而矩阵A是实对称矩阵，所有特征值均为实数，故我们可以用幂法和反幂法结合适当的平移量求解。其中求解最大的特征值和最小的特征值，采用幂法并结合适当的平移量求解，而求解λs，和距离μk最近的特征值采用带平移的反幂法求解。算法分析一、幂法求解λ我们知道幂法适用于求解矩阵A的绝对值最大的特征值，设矩阵A的特征值为λ1≥λ2≥⋯≥|λn|，对于情况λ1>λ2，幂法是适用的，并可以求出特征值λ1，然而对于λ此时，由于λ1≤λ2≤⋯≤λ501B=A-λI然后对矩阵B用幂法。由线性代数知识知，矩阵B的特征值为λ1-λ,λ2-λ,⋯,λ501-λ，若λ=λ1，则有0=λ1-λ≤λ2-λ≤⋯≤λ501-λ，对B综上，对于λ1≠-λ501，只需要对A用幂法即可得到A的绝对值最大的特征值λ，然后做平移B=A-λI，对B用幂法即可求λ由于这种情况，不能直接对矩阵A使用幂法。此时有两种解决办法，一种是做适当平移B=A-pI，对用幂法迭代求出其绝对值最大的特征值λ后，则-P+λ，P+λ={λ综合1，2，求解矩阵A的特征值λ1,λ501算法思路为设定最大迭代次数和收敛精度，然后对矩阵A使用幂法，如果收敛，说明λ1≠-λ501，则可以得到一个特征值λ，然后做平移B=A-λI，对矩阵B使用幂法得到B的一个特征值λ'，则有{λ1,λ501}=λ,λ'+λ；若达到最大迭代次数还没有达到收敛要求的精度，说明λ二、反幂法求解λ设Bμ=A-μI-1，由于Bμ的特征值为1λ1-μ,1λ2-μ,⋯,1λ501-μ，可知以次取μ=0,u1,μ2,⋯,μ39时，矩阵Bμ三、求cond由定义知，condA2=A2A-12=λmλs,λm=max1≤i≤501λi，λm,λs在前面都已经求出，可以直接用来算出condA2。而对于最后对于矩阵A的存储，考虑到矩阵A的特殊性，只需定义一个一维数组用来存储矩阵A的对角线上的元素，和两个浮点数存储b，c即可。算法实现流程图幂法求解λ1反幂法求λs计算结果λ-1.070011361502E+01λ-9.935586060075E-01λ9.724634099619E+00λ-4.870426738850E-01λ-5.557910794229E-03λ2.231736249575E-02λ-1.018293403315E+01λ5.324174742069E-01λ-9.585707425068E+00λ1.052898962693E+00λ-9.172672423928E+00λ1.589445881881E+00λ-8.652284007898E+00λ2.060330460274E+00λ-8.093483808675E+00λ2.558075597073E+00λ-7.659405407692E+00λ3.080240509307E+00λ-7.119684648691E+00λ3.613620867692E+00λ-6.611764339397E+00λ4.091378510451E+00λ-6.066103226595E+00λ4.603035378279E+00λ-5.585101052628E+00λ5.132924283898E+00λ-5.114083529812E+00λ5.594906348083E+00λ-4.578872176865E+00λ6.080933857027E+00λ-4.096470926260E+00λ6.680354092112E+00λ-3.554211215751E+00λ7.293877448127E+00λ-3.041090018133E+00λ7.717111714236E+00λ-2.533970311130E+00λ8.225220014050E+00λ-2.003230769563E+00λ8.648666065193E+00λ-1.503557611227E+00λ9.254200344575E+00cond1.925204273902e+003det2.772786141752e+118初始迭代向量为：u0=(1,1,1,⋯,1)。分析初始向量对结果的影响分析初始向量对计算结果的影响，以计算λ1,λ501为例(初始迭代向量λ1迭代次数正确与否是否收敛u-1.070011361502e+001341是是u-2.080981085336e+000158否是μ-2.080981085338e+000184否是u-1.070011361509e+001412是是初始迭代向量λ50迭代次数正确与否是否收敛μ9.723340141234e+00010000是否u9.724634098771e+000641是是u9.724634099653e+00081否是u9.724634099652e+00073是是u9.724634098776e+0001023是是分析上表的结果可以看到，使用幂法求矩阵特征值的时候，初始向量的选取对结果的影响：一方面，会影响幂法的收敛速度，使得迭代收敛缓慢，如取初始向量为u0=(1,1,1,⋯,1)计算λ501时，迭代次数达到给定的最大值时仍未收敛；另一方面可能会得不到想要的结果，如取u0程序源代码幂法#include<stdio.h>#include<math.h>#include<stdlib.h>doubleerro(doublelamd0,doublelamd1);//每步迭代结束后计算误差的函数voiddiedai(doublea[],doubleb,doublec,double*x0,double*x1);//幂法迭代函数doubledianji(doublex[],doubley[]);//计算两个向量内积的函数intmain(){ doublea[501]={0},x[505]={0.0},y[505]={0.0},b=0.16,c=-0.064,lamd=0,lamd0=0,lamd1=0,e1=0,d=0,*x1=y,*x0=x,*t; inti,k1,k2;//变量定义即初始化，其中一维数组x，y用于存储每步迭代后的向量，指针变量x0，x1用于存取x，y的地址用于迭代计算。为了迭代格式的统一，对初始向量首尾各扩充了两个零元素 for(inti=1;i<502;++i) { a[i-1]=(1.64-i*0.024)*sin(0.2*i)-0.64*exp(0.1/i); x[i+1]=1;//赋值初始迭代向量 } for(i=0;i<10000;++i) { d=sqrt(dianji(x0,x0));//计算初始向量也即uk-1的2范数 for(intj=2;j<503;j++) { x0[j]=x0[j]/d;//计算yk-1 } diedai(a,b,c,x0,x1);//根据向量yk-1计算x1，uk lamd1=dianji(x0,x1);//计算yk-1，uk的内积，即βk e1=erro(lamd0,lamd1);//计算误差 if(e1<1e-12)//判断是否收敛 { lamd=lamd1; break; } else//未收敛则将uk，βk用于下次迭代 t=x0; x0=x1; x1=t; lamd0=lamd1; } k1=i-1; for(inti=0;i<501;++i) { a[i]=a[i]-lamd;//以第一次计算结果为平移量进行幂法迭代 x[501+1]=i; } for(i=0;i<100000;++i) { d=sqrt(dianji(x0,x0)); for(intj=2;j<503;j++) { x0[j]=x0[j]/d; } diedai(a,b,c,x0,x1); lamd1=dianji(x0,x1); e1=erro(lamd0,lamd1); if(e1<1e-12) { lamd1=lamd1+lamd; k2=i;//记录迭代步数 break; } else t=x0; x0=x1; x1=t; lamd0=lamd1; } lamd1=lamd1+lamd; k2=i; printf("lamd1=%20.12e,K1=%d\nlamd501=%20.12e,K2=%d\n",lamd<lamd1?lamd:lamd1,lamd<lamd1?k1:k2,lamd>lamd1?lamd:lamd1,lamd>lamd1?k1:k2);//输出结果 system("pause"); return0;}doubledianji(doublex[],doubley[])//定义向量内积{ doublee=0; for(inti=2;i<503;++i) { e=e+x[i]*y[i]; } return(e);}doubleerro(doublelamd0,doublelamd1)//计算两个数之间的相对误差{ doublee,d,l; e=fabs(lamd1-lamd0); d=fabs(lamd1); l=e/d; return(l);}voiddiedai(doublea[],doubleb,doublec,double*x0,double*x1){ for(inti=0;i<501;++i) { x1[i+2]=c*x0[i]+b*x0[i+1]+a[i]*x0[i+2]+b*x0[i+3]+c*x0[i+4];//计算uk=Ayk-1 }}反幂法#include<stdio.h>#include<math.h>#include<stdlib.h>doubleerro(doublelamd0,doublelamd1);//每步迭代结束后计算误差的函数voiddiedai(double(*L)[501],double(*U)[3],double*x0,double*x1);//反幂法迭代函数，根据LU分解函数结果反解ukdoubledianji(doublex[],doubley[]);//计算两个向量内积的函数voidLU(doublea[],doubleb,doublec,double(*L)[501],double(*U)[3]);//对A进行LU分解intmain(){ doublea[501]={0},aa[501],x[501]={0.0},y[501]={0.0},l[3][501]={0},u[501][3]={0},(*L)[501],(*U)[3],b=0.16,c=-0.064,lamd0=0,lamd1=0,e1=0,d=0,*x1=y,*x0=x,*t; for(inti=1;i<502;++i)//变量定义即初始化，其中一维数组x，y用于存储每步迭代后的向量，指针变量x0，x1用于存取x，y的地址用于迭代计算。其中3X501维数组l,u用于存放对a分解后的矩阵L，U，指针数组L，U用于存放l，u的元素地址 { a[i-1]=(1.64-i*0.024)*sin(0.2*i)-0.64*exp(0.1/i); x[i-1]=1;//赋值初始迭代向量 } L=l; U=u; x0=x; x1=y;LU(a,b,c,L,U);//对A进行LU分解 for(inti=0;i<1000;++i) { d=sqrt(dianji(x0,x0));//计算初始向量也即uk-1的2范数 for(intj=0;j<501;j++) { x0[j]=x0[j]/d;//计算yk-1 } diedai(L,U,x0,x1);/根据向量yk-1反解x1即uk lamd1=dianji(x0,x1);//计算yk-1，uk的内积，即βk e1=erro(lamd0,lamd1);/计算误差 if(e1<1e-12)//判断是否收敛 { printf("lamds=%.12e,K=%d\n",1/lamd1,i); break; } else t=x0;//未收敛则将uk，βk用于下次迭代 x0=x1; x1=t; lamd0=lamd1; } doubledetA=1,condA,lamds=1/lamd1; for(inti=0;i<501;++i) { detA=detA*U[i][0];//根据LU分解所得的矩阵U计算detA } printf("detA=%.12e\n",detA); doublelamd=-10.700113615016,lamd501=9.724634099619,delta,miu; condA=fabs(lamd/lamds); printf("condA=%.12e\n",condA); delta=(lamd501-lamd)/40; for(intk=1;k<40;++k)//计算距离μk最近的特征值 { miu=lamd+delta*k;//计算μk for(inti=0;i<501;++i) { aa[i]=a[i]-miu;//对矩阵A进行平移量为μk的平移 x[i]=1; y[i]=2; } LU(aa,b,c,L,U);//计算平移后的矩阵A的LU分解，并存储到数组l，u中用于下面的迭代计算 lamd0=0; for(inti=0;i<1000;++i) { d=sqrt(dianji(x0,x0)); for(intj=0;j<501;j++) { x0[j]=x0[j]/d; } diedai(L,U,x0,x1); lamd1=dianji(x0,x1); e1=erro(lamd0,lamd1); if(e1<1e-12) { printf("lamdi_%d=%.12e,K=%d\n",k,1/lamd1+miu,i); break; } else t=x0; x0=x1; x1=t; lamd0=lamd1; } } system("pause"); return0;}doubledianji(doublex[],doubley[])//定义向量内积{ doublee=0; for(inti=0;i<501;++i) { e=e+x[i]*y[i]; } return(e);}doubleerro(doublelamd0,doublelamd1){ doublee,d,l; e=fabs(lamd1-lamd0); d=fabs(lamd1); l=e/d; return(l);}voidLU(doublea[],doubleb,doublec,double(*L)[501],double(*U)[3])//计算矩阵A类型的矩阵的LU分解{ U[0][0]=a[0]; U[0][1]=b; U[0][2]=c; L[1][0]=b/U[0][0]; L[2][0]=c/U[0][0]; U[1][0]=a[1]-L[1][0]*U[0][1]; U[1][1]=b-L[1][0]*U[0][2]; U[1][2]=c;L[1][1]=(b-L[2][0]*U[0][1])/U[1][0];L[2][1]=c/U[1][0]; for(intk=2;k<499;++k) { U[k][0]=a[k]-L[2][k-2]*U[k-2][2]-L[1][k-1]*U[k-1][1]; U[k][1]