北航数值分析大作业一.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矩阵中的零元素，只保存主对角元素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> double ep=1e-12,b=0.16,c=-0.064; int j=0; double power(double a[501]); //幂法 double inv_power(double a[501]); //反幂法 double det(double a[501]) ; //求det int main() //主程序 { int i,k; double A[501],B[501],beta_1,beta_501,beta_s,beta_k; double mu; 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)); } double power(double a[501]) //幂法 { int i=0,N=5000; double b=0.16,c=-0.064; double u[501],y[501]; double m=1,beta;