北航数值分析大作业一

北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号ZY1403140学生姓名许阳教师孙玉泉日期2014年11月26日设有501501的实对称矩阵A，5011Aabcbccbcba其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0cbieiiaii。矩阵A的特征值为)501,,2,1(ii，并且有||min||,501150121iis1.求1，501和s的值。2.求A的与数4015011kk最接近的特征值)39,,2,1(kki。3.求A的(谱范数)条件数2)A(cond和行列式detA。一方案设计1求1，501和s的值。s为按模最小特征值，||min||5011iis。可使用反幂法求得。1，501分别为最大特征值及最小特征值。可使用幂法求出按模最大特征值，如结果为正，即为501，结果为负，则为1。使用位移的方式求得另一特征值即可。2求A的与数4015011kk最接近的特征值)39,...,2,1(kki。题目可看成求以k为偏移量后，按模最小的特征值。即以k为偏移量做位移，使用反幂法求出按模最小特征值后，加上k，即为所求。3求A的(谱范数)条件数2)(Acond和行列式detA。矩阵A为非奇异对称矩阵，可知，||)(minmax2Acond(1-1)其中max为按模最大特征值，min为按模最小特征值。detA可由LU分解得到。因LU均为三角阵，则其主对角线乘积即为A的行列式。二算法实现1幂法使用如下迭代格式：||max|)|sgn(max||max/),,(111111)0()0(10kkkkkkkkTnuuAyuuuyuuu任取非零向量(2-1)终止迭代的控制理论使用||/||1kkk，实际使用||/||||||1kkk(2-2)由于不保存A矩阵中的零元素，只保存主对角元素a[501]及b,c值。则上式中1kkAyu简化为：)499,,3()2()1()()()1()2()()501()501()500()499()501()501()500()500()499()498()500()4()3()2()2()1()2()3()2()1()1()1(iicyibyiyiaibyicyiuyabycyubyyabycyucybyyabyucybyyau(2-3)2反幂法使用如下迭代格式：||max|)|sgn(max||max/),,(1111-111)0()0(10kkkkkkkkTnuuyAuuuyhhu任取非零向量(2-4)其中111kkkkyAuyAu，解方程求出ku。求解过程中使用LU分解，由于A为5对角矩阵，选择追赶法求取LU分解。求解过程如下：kkkkkkkuxUuyLxyLUu11追赶法求LU分解的实现：11500499115015015013215011tqqtprzzrpLUabcbccbcbaA(2-5)由上式推出分解公式如下：501,...,3,501,...,3,501,...,3,500,...,2,/)(499,...,1,/,/,1221122221111itrcqapictbriczipqrbtipcqtrapbrabtapiiiiiiiiiiiiii(2-6)推导出回代求解公式如下：501,...,3,/)(/)(/1221222111ipxrxzyxpxryxpyxiiiiiii(2-7)1,...,499,21501500500500501501ixqutxuutxuxuiiiiii(2-8)32)A(cond及A行列式求解||||)A(12scond(2-9)由式(2-5)可得：5011detiipA(2-10)三源程序#includestdio.h#includemath.hdoubleep=1e-12,b=0.16,c=-0.064;intj=0;doublepower(doublea[501]);//幂法doubleinv_power(doublea[501]);//反幂法doubledet(doublea[501]);//求detintmain()//主程序{inti,k;doubleA[501],B[501],beta_1,beta_501,beta_s,beta_k;doublemu;for(i=0;i501;i++)A[i]=(1.64-0.024*(i+1))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1));beta_1=power(A);//第一问printf(λ1\t=%.12e\t迭代次数：%d\n,beta_1,j);for(i=0;i501;i++)//位移B[i]=A[i]-beta_1;beta_501=power(B)+beta_1;printf(λ501\t=%.12e\t迭代次数：%d\n,beta_501,j);beta_s=inv_power(A);printf(λs\t=%.12e\t迭代次数：%d\n,beta_s,j);for(k=1;k=39;k++)//第二问{mu=beta_1+k*(beta_501-beta_1)/40;for(i=0;i501;i++)B[i]=A[i]-mu;beta_k=inv_power(B)+mu;printf(λi%d\t=%.12e\t迭代次数：%d\n,k,beta_k,j);}printf(cond(A)2=%.12e\n,beta_1/beta_s);//第三问printf(detA\t=%.12e\n,det(A));}doublepower(doublea[501])//幂法{inti=0,N=5000;doubleb=0.16,c=-0.064;doubleu[501],y[501];doublem=1,beta;for(i=0;i501;i++)u[i]=1;j=0;while(jN){for(i=0;i501;i++){y[i]=u[i]/fabs(m);}u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3];u[499]=c*y[497]+b*y[498]+a[499]*y[499]+b*y[500];u[500]=c*y[498]+b*y[499]+a[500]*y[500];for(i=2;i499;i++){u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}beta=0;for(i=0;i501;i++){if(fabs(u[i])=fabs(beta))beta=u[i];}if(beta0)if(fabs(fabs(beta)-fabs(m))/fabs(beta)ep)break;if(fabs(beta-m)/fabs(beta)ep)break;m=beta;j++;}returnbeta;}doubleinv_power(doublea[501])//反幂法{doublep[501],r[501],t[501],q[501],u[501],y[501];doublebeta,m=1;inti,N=1000;p[0]=a[0];t[0]=b/p[0];r[1]=b;p[1]=a[1]-r[1]*t[0];q[0]=c/p[0];q[1]=c/p[1];t[1]=(b-r[1]*q[0])/p[1];for(i=2;i501;i++){r[i]=b-c*t[i-2];p[i]=a[i]-c*q[i-2]-r[i]*t[i-1];q[i]=c/p[i];t[i]=(b-r[i]*q[i-1])/p[i];}for(i=0;i501;i++)u[i]=1;j=0;while(jN){for(i=0;i501;i++){y[i]=u[i]/fabs(m);}u[0]=y[0]/p[0];u[1]=(y[1]-r[1]*u[0])/p[1];for(i=2;i501;i++)u[i]=(y[i]-c*u[i-2]-r[i]*u[i-1])/p[i];u[499]=u[499]-t[499]*u[500];for(i=498;i=0;i--)u[i]=u[i]-t[i]*u[i+1]-q[i]*u[i+2];beta=0;for(i=0;i501;i++){if(fabs(u[i])=fabs(beta))beta=u[i];}if(beta0)if(fabs(fabs(beta)-fabs(m))/fabs(beta)ep)break;if(fabs(beta-m)/fabs(beta)ep)break;m=beta;j++;}return1/beta;}doubledet(doublea[501])//求det{doubledet_A=1;doublep[501],r[501],t[501],q[501];inti;p[0]=a[0];t[0]=b/p[0];r[1]=b;p[1]=a[1]-r[1]*t[0];q[0]=c/p[0];q[1]=c/p[1];t[1]=(b-r[1]*q[0])/p[1];for(i=2;i501;i++){r[i]=b-c*t[i-2];p[i]=a[i]-c*q[i-2]-r[i]*t[i-1];q[i]=c/p[i];t[i]=(b-r[i]*q[i-1])/p[i];}for(i=0;i501;i++)det_A=det_A*p[i];returndet_A;}四程序结果五计算过程中的现象使用||/||1kkk作为终止迭代条件时，出现迭代无法终止的情况，通过调试发现按模最大特征值为负时，当k充分大后，迭代向量ku各分量不断变号，使得k与1k异号，判别式||/||1kkk不收敛。因此将终止迭代条件修改为||/||||||1kkk，程序实现如下：if(beta0)if(fabs(fabs(beta)-fabs(m))/fabs(beta)ep)break;if(fabs(beta-m)/fabs(beta)ep)break;从迭代次数可以看出1与501收敛较慢，由按模最大特征值与按模次大特征值的比值越小，收敛速度越慢，可知存在与1和501的模相近的特征值。