北航数值分析大作业一.docx
北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号ZY1403140学生姓名许阳教师孙玉泉日期2021年11月26日设？!501x501的实对称矩阵A,0.1其中，ai=(1.64-0.024z)sin(0.2z)-0.64e(i=1,2,501),Z?=0.16,c=-0.064。矩阵A的特征值为4(i=l,2,.,501),并且有4，4o和儿的值。A的与数4=4+A4。；4最接近的特征值ik(k=1,2,39)。A的(谱范数)条件数CaM（八）2和行列式CletAo一方案设计1求4，4(H和4的值。4为按模最小特征值，"=min141。可使用反幕法求得。l50l1,Zoi分别为最大特征值及最小特征值。可使用事法求出按模最大特征值，如结果为正，即为4结果为负，那么为乙。使用位移的方式求得另一特征值即可。2求A的与数4=4+14o-4最接近的特征值i(k=1,2,.,39)。40题目可看成求以4为偏移量后，按模最小的特征值。即以4为偏移量做位移，使用反吊法求出按模最小特征值后，加上以，即为所求。3求A的(谱范数)条件数CMd（八）2和行列式detAo矩阵A为非奇异对称矩阵，可知,Cond(A)2 =|4nax4nin(1-1)其中4”为按模最大特征值，AnM为按模最小特征值。detA可由LlJ分解得到。因LU均为三角阵，那么其主对角线乘积即为A的行列式。二算法实现1幕法使用如下迭代格式：任取非零向勤O=(Y,4。Y加=kmaXlkl(21)uk=A%k=sgn(maxw-11)max|以Tl终止迭代的控制理论使用I瓦-片I/1凤£,实际使用IIAI-IA-Ivk<(2-2)由于不保存A矩阵中的零元素，只保存主对角元素a501及b,c值。那么上式中uk=Ayj简化为:M(I)=a(l)y(l)+y(2)+Cy(3)u(T)=by(J)+(2)y(2)+by(3)+Cy(4)w(500)=cy(498)+by(499)+(500)y(500)+切(501)”(501)=Cy(499)+力(500)+4(501)y(501)u(i)=cy(i-2)÷by(i-1)+()y()+by(i+1)+cy(i+2)(i=3,499)2反暴法使用如下迭代格式：任取非零向羸O=(AI(°),冏°)尸%=以/maxI/Tl(24)uk=AyTk=sgn(maxuk_1)maxuk_x其中uk=Aiyk-i=>Auk=%_,解方程求出以O求解过程中使用LU分解，由于A为5对角矩阵，选择追赶法求取LU分解。求解过程如下:追赶法求LU分解的实现：ahPl r2=Z3由上式推出分解公式如下：(2-5)A=c推导出回代求解公式如下:PI/=b!ax弓="p2=a2-r2tl¾=cpf.,f=l,.499Ti=S-d-)Pi,i=2,500(2-6)zj=c,z=3,.5lri=/?-cZf_2,f=3,.501Pi=ai-Cq_2-=3,.501玉二yP<2=(J2-¾)P2(2-7)=(M_zix-2-)ppi=3,.501m50I=毛01"m5=x500-500m501(2-8)%=XiTMH-仇+2，i=499,J3COnd（八）2及A行列式求解由式(2-5)可得:三源程序#include<stdio.h>#include<math.h>cond（八）2=-501detA=11pi1=1(2-9)(2-10)doubleep=le-12,b=0.16,c=-0.064;intj=O;doublepower(doublea501);嘉法doubleinv_power(doublea501);反事法doubledet(doublea501);求detint main()主程序inti,k;doubleA501,B501,beta_1,beta_501,beta_s,beta_k;doublemu;for(i=0;i<501;i+)Ai=(1.64-0.024*(i+l)*sin(0.2*(i+l)-0.64*exp(0.1(i+1);beta_1=power(A);第一问printf(',lt=%.12et迭代次数：%dn",beta.l,j);for(i=0;i<501;i+)位移Bli=Ali-beta_l;beta_501=power(B)+beta_l;printf<',501t=%.12et迭代次数：%dn',beta.501,j);beta_s=inv_power（八）;PrintfCst=%12et迭代次数：%dn",beta_s,j);for(k=1;k<=39;k+)第二问(mu=beta_l+k*(beta-501-beta_1)/40;for(i=0;i<501;i+)Bi=Ali-mu;beta_k=inv_power(B)+mu;printf(,i%dt=%.12et迭代次数：%dnn,k,beta_k,j);)printf(,cond（八）2=%.12en",beta_l/beta_s);/第三问printf(',detAt=%.12en",det（八）);)doublepower(doublea501)事法(inti=0,N=5000;doubleb=0.16,c=-0.064;doubleul501,y501;doublem=l,beta;for(i=0;i<501;i+)ui=l;j=0;while(j<N)(for(i=0;i<501;i+)yli=uilfabs(m);uO=aO*yO+b*y1+c*yl2;u1=b*yO+a1*y1+b*y2+c*y3;u499=c*y497+b*y498+a4991*y499+b*yl500;u500=c*yl498+b*y499+a500*y500;for(i=2;i<499;i+)ui=c*yi-2+b*yli-l+ai*yi+b*yi+1+c*yi+2;beta=O;for(i=0;i<501;i+)if(fabs(ui)>=fabs(beta)beta=ui;if(beta<O)if(fabs(fabs(beta)-fabs(m)fabs(beta)<ep)break;if(fabs(beta-m)fabs(beta)<ep)break;m=beta;j+;)returnbeta;)doubleinv_power(doublea5Ol)反嘉法(doublep501,r501,t501,q501,u501,y501;doublebeta,m=l;inti,N=100O;pO=aO;tO=b/pO;rl=b;pl=al-rl*tO;qO=c/pO;ql=c/pl;tl=(b-rl*qO)pl;for(i=2;i<501;i+)(rli=b-c*tli-2J;pi=ai-c*qli-2-ri*ti-1;qi=cpi;ti=(b-ri*qi-l)pi;)for(i=0;i<501;i+)ui=l;j=0；while(j<N)for(i=0;i<501;i+)yi=uifabs(m);uO=yOpO;ul=(yl-rl*uO)pl;for(i=2;i<501;i+)ui=(yi-c*ui-2-ri*ui-1)pi;u499=u499-t499*u500;for(i=498;i>=0;i-)uil=ui-ti*ui+l-qli*ui+2;beta=O;for(i=0;i<501;i+)(if(fabs(ui)>=fabs(beta)beta=ui;)if(beta<O)if(fabs(fabs(beta)-fabs(m)fabs(beta)<ep)break;if(fabs(beta-m)fabs(beta)<ep)break;m=beta;j+;)return1/beta;)doubledet(doublea501)求detdoubledet_A=1;doublep50l,r501,t501,q501;inti;pO=aO;tO=b/pO;rl=b;pl=al-rl*tO;q0=c/p0;ql=c/pl;tU=(b-rl*qO)pl;for(i=2;i<501;i+)(ri=b-c*ti-2;pi=ai-c*qi-2-ri*ti-1;qi=cpi;ti=(b-ri*qi-l)pi;)for(i=0;i<501;i+)det_A=det_A*pi;returndet_A;)四程序结果五计算过程中的现象使用1瓦-瓦T|/|凤lg作为终止迭代条件时，出现迭代无法终止的情况，通过调试发现按模最大特征值为负时，当k充分大后，迭代向量人各分量不断变号，使得用与AT异号，判别式I瓦-ATI/IAl不收敛。因此将终止迭代条件修改为IlAI-IA-.Il/1A£，程序实现如下：if(beta<O)if(fabs(fabs(beta)-fabs(m)fabs(beta)<ep)break;if(fabs(beta-m)fabs(beta)<ep)break;从迭代次数可以看出4与4OI收敛较慢，由按模最大特征值与按模次大特征值的比值越小，收敛速度越慢，可知存在与4和4(H的模相近的特征值。【本文档内容可以自由复制内容或自由编辑修改内容期待你的好评和关注，我们将会做得更好】
