北航数值分析大作业一.doc

北 京 航 空 航 天 大 学数值分析大作业一学院名称 自动化 专业方向 控制工程 学 号 ZY1403140 学生姓名 许阳 教师 孙玉泉 日 期 2014 年 11月26 日 设有的实对称矩阵A，其中，。矩阵A的特征值为，并且有1.求，和的值。2.求A的与数最接近的特征值。3.求A的(谱范数)条件数和行列式detA。一 方案设计1 求，和的值。 为按模最小特征值，。可使用反幂法求得。 ，分别为最大特征值及最小特征值。可使用幂法求出按模最大特征值，如结果为正，即为，结果为负，则为。使用位移的方式求得另一特征值即可。2 求A的与数最接近的特征值。 题目可看成求以为偏移量后，按模最小的特征值。即以为偏移量做位移，使用反幂法求出按模最小特征值后，加上，即为所求。3 求A的(谱范数)条件数和行列式detA。 矩阵A为非奇异对称矩阵，可知，(1-1) 其中为按模最大特征值，为按模最小特征值。 detA可由LU分解得到。因LU均为三角阵，则其主对角线乘积即为A的行列式。二 算法实现1 幂法 使用如下迭代格式：(2-1) 终止迭代的控制理论使用， 实际使用(2-2) 由于不保存A矩阵中的零元素，只保存主对角元素a501及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 #include double ep=1e-12,b=0.16,c=-0.064; int j=0;double power(double a501); /幂法double inv_power(double a501); /反幂法 double det(double a501) ; /求det int main()/主程序int i,k;double A501,B501,beta_1,beta_501,beta_s,beta_k;double mu;for(i=0;i501;i+)Ai=(1.64-0.024*(i+1)*sin(0.2*(i+1)-0.64*exp(0.1/(i+1); beta_1=power(A); /第一问printf(1t= %.12et迭代次数：%dn,beta_1,j);for(i=0;i501;i+) /位移 Bi=Ai-beta_1; beta_501=power(B)+beta_1; printf(501t= %.12et迭代次数：%dn,beta_501,j); beta_s=inv_power(A); printf(st= %.12et迭代次数：%dn,beta_s,j); for(k=1;k=39;k+) /第二问 mu=beta_1+k*(beta_501-beta_1)/40; for(i=0;i501;i+) Bi=Ai-mu; beta_k=inv_power(B)+mu; printf(i%dt= %.12et迭代次数：%dn,k,beta_k,j); printf(cond(A)2= %.12en,beta_1/beta_s); /第三问printf(detAt= %.12en,det(A);double power(double a501) /幂法 int i=0,N=5000;double b=0.16,c=-0.064;double u501,y501;double m=1,beta;for(i=0;i501;i+) ui=1; j=0;while(jN)for(i=0;i501;i+)yi=ui/fabs(m); u0=a0*y0+b*y1+c*y2; u1=b*y0+a1*y1+b*y2+c*y3; u499=c*y497+b*y498+a499*y499+b*y500; u500=c*y498+b*y499+a500*y500; for(i=2;i499;i+) ui=c*yi-2+b*yi-1+ai*yi+b*yi+1+c*yi+2; beta=0; for(i=0;i=fabs(beta) beta=ui; if(beta0) if(fabs(fabs(beta)-fabs(m)/fabs(beta)ep) break; if(fabs(beta-m)/fabs(beta)ep) break; m=beta;j+; return beta;double inv_power(double a501) /反幂法 double p501,r501,t501,q501,u501,y501;double beta,m=1;int i,N=1000;p0=a0;t0=b/p0;r1=b;p1=a1-r1*t0;q0=c/p0;q1=c/p1;t1=(b-r1*q0)/p1;for(i=2;i501;i+)ri=b-c*ti-2;pi=ai-c*qi-2-ri*ti-1;qi=c/pi;ti=(b-ri*qi-1)/pi;for(i=0;i501;i+)ui=1;j=0;while(jN)for(i=0;i501;i+)yi=ui/fabs(m); u0=y0/p0;u1=(y1-r1*u0)/p1;for(i=2;i=0;i-)ui=ui-ti*ui+1-qi*ui+2;beta=0; for(i=0;i=fabs(beta) beta=ui; if(beta0) if(fabs(fabs(beta)-fabs(m)/fabs(beta)ep) break; if(fabs(beta-m)/fabs(beta)ep) break; m=beta;j+; return 1/beta;double det(double a501) /求det double det_A=1;double p501,r501,t501,q501;int i;p0=a0;t0=b/p0;r1=b;p1=a1-r1*t0;q0=c/p0;q1=c/p1;t1=(b-r1*q0)/p1;for(i=2;i501;i+)ri=b-c*ti-2;pi=ai-c*qi-2-ri*ti-1;qi=c/pi;ti=(b-ri*qi-1)/pi;for(i=0;i501;i+)det_A=det_A*pi;return det_A;四 程序结果 五 计算过程中的现象 使用 作为终止迭代条件时，出现迭代无法终止的情况，通过调试发现按模最大特征值为负时，当k充分大后，迭代向量各分量不断变号，使得与异号，判别式不收敛。因此将终止迭代条件修改为，程序实现如下：if