北航数值分析大作业一

北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号ZY1403140学生姓名许阳教 师孙玉泉日期2021年11月26日设有的实对称矩阵A，其中，。矩阵A的特征值为，并且有，和的值。A的与数最接近的特征值。A的(谱范数)条件数和行列式detA。一方案设计1求，和的值。为按模最小特征值，。可使用反幂法求得。，分别为最大特征值及最小特征值。可使用幂法求出按模最大特征值，如结果为正，即为，结果为负，那么为。使用位移的方式求得另一特征值即可。2求A的与数最接近的特征值。题目可看成求以为偏移量后，按模最小的特征值。即以为偏移量做位移，使用反幂法求出按模最小特征值后，加上，即为所求。3求A的(谱范数)条件数和行列式detA。矩阵A为非奇异对称矩阵，可知， (1-1)其中为按模最大特征值，为按模最小特征值。detA可由LU分解得到。因LU均为三角阵，那么其主对角线乘积即为A的行列式。二算法实现1幂法使用如下迭代格式： (2-1)终止迭代的控制理论使用，实际使用 (2-2)由于不保存A矩阵中的零元素，只保存主对角元素a[501]及b,c值。那么上式中简化为： (2-3)2反幂法使用如下迭代格式： (2-4)其中，解方程求出。求解过程中使用LU分解，由于A为5对角矩阵，选择追赶法求取LU分解。求解过程如下：追赶法求LU分解的实现： (2-5)由上式推出分解公式如下： (2-6)推导出回代求解公式如下： (2-7) (2-8)3及A行列式求解 (2-9)由式(2-5)可得： (2-10)三源程序#include<stdio.h>#include<math.h>doubleep=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;i<501;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;i<501;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;i<501;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;i<501;i++) u[i]=1;j=0; while(j<N) { for(i=0;i<501;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;i<499;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;i<501;i++) { if(fabs(u[i])>=fabs(beta)) beta=u[i]; } if(beta<0) 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;i<501;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;i<501;i++) u[i]=1; j=0; while(j<N) { for(i=0;i<501;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;i<501;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;i<501;i++) { if(fabs(u[i])>=fabs(beta)) beta=u[i]; } if(beta<0) 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;i<501;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;i<501;i++) det_A=det_A*p[i]; returndet_A;}四程序结果五计算过程中的现象使用作为终止迭代条件时，出现迭代无法终止的情况，通过调试发现按模最大特征值为负时，当k充分大后，迭代向量各分量不断变号，使得与异号，判别式不收敛。因此将终止迭代条件修改为，程序实现如下： if(beta<0)