北航数值分析大作业一

1、数值分析第一次作业设有的实对称矩阵A，其中，。矩阵A的特征值为，并且有1.求，和的值。2.求A的与数最接近的特征值。3.求A的(谱范数)条件数和行列式detA。一 方案设计1 求，和的值。 是矩阵A最小的特征值，为矩阵最大的特征值。对矩阵A使用幂法可以求出、其中的一个，但不确定是哪一个，设求出的值为。然后再构造矩阵B=A-I，对矩阵B使用幂法求出其按模最大的特征值，可得、中另外一个值，为+。比较和+，其中较小者为为，较大者为。为按模最小特征值，,直接利用反幂法可以求得。2 求A的与数最接近的特征值。 题目可看成求以为偏移量后，按模最小的特征值。即以为偏移量做位移，使用反幂法求出按模最小特征值后

2、，加上，即为所求。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分解法和JACOBI迭代法两种。 （1）LU分解由于A为5对角矩阵，选择追赶法求取LU分解。求解过程如

3、下：(2-5) 由上式推出分解公式如下：(2-6)推导出回代求解公式如下：(2-7)(2-8)（2）JACOBI迭代法迭代公式如下所示：u1k=-bu2k-1+cu3k-1+y1k-1/a11u2k=-(bu1k-1+bu3k-1cu4k-1)+ y2k-1/a22u500k=-bu498k-1+bu499k-1+cu501k-1+y500k-1/a500500 u501k=-cu499k-1+bu500k-1+y501k-1/a501501uik=-cui-2k-1+bui-1k-1+bui+1k-1+cui+2k-1+yik-1/aii(2-9)3 及A行列式求解 (2-10) 由式(2-

4、5)可得： (2-11)三 源程序#include #include double mifa(double a501); /幂法double fanmifa1(double a501); /反幂法（三角分解法求解U(k)） double det(double a501) ; /求det double fanmifa(double a501);/反幂法（jacobi求解）double ep=1e-12,b=0.16,c=-0.064; int j=0;int main()int i,k;double A501,B501,beta_1,beta_501,beta_s,beta_k,beta_m;d

5、ouble 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=mifa(A); for(i=0;ibeta_1)beta_1=beta_1;beta_501=beta_m;elsebeta_1=beta_m;beta_501=beta_1;printf(1t= %.12etn,beta_1);for(i=0;i501;i+) Bi=Ai-beta_1; beta_501=mifa(B)+beta_1; printf(501=t= %.12etn,beta_501); beta_s=f

6、anmifa(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=fanmifa(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);return 0;double mifa(double a501) /幂法 int i

7、=0,r;double b=0.16,c=-0.064;double u501,y501;double s,lastbeta,beta,sq;for(i=0;i501;i+) ui=1;j=0;while(1)j+;s=fabs(u0);for(i=1;i501;i+)sq=fabs(ui);if(ssq) s=sq;for(i=0;i501;i+)if(fabs(ui)=s) r=i;for(i=0;i501;i+)yi=ui/s;u0=a0*y0+b*y1+c*y2; /迭代求出下一个u u1=b*y0+a1*y1+b*y2+c*y3; u499=c*y497+b*y498+a499*y4

8、99+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;lastbeta=beta;for(i=0;i501;i+)beta=ur/yr; if(fabs(beta-lastbeta)/fabs(beta)ep) break; return beta;double fanmifa1(double a501) /反幂法(用三角分解法求解AU(k)=Y(k-1)中的U(k) double p501,r501,t501,q501,u501,y501;double bet

9、a,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

10、;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 fanmifa(double a501) /用jacobi迭代法求解的反幂法double u501,y501,lastu501,beta=0,lastbeta;double s=0,max;int i,r;for(i=0;i501;i+) /将初始向量的分量全赋值为1ui=1;j=0;while(1

11、)j+;s=0;for(i=0;i501;i+)s=ui*ui+s;s=sqrt(s);for(i=0;i501;i+) /将y单位化yi=ui/s;while(1)for(i=0;i501;i+)lastui=ui;u0=(-b*lastu1-c*lastu2+y0)/a0;u1=(-b*lastu0-b*lastu2-c*lastu3+y1)/a1;u499=(-c*lastu497-b*lastu498-b*lastu500+y499)/a499;u500=(-c*lastu498-b*lastu499+y500)/a500;for(i=2;i499;i+)ui=(-c*lastui-2

12、-b*lastui-1-b*lastui+1-c*lastui+2+yi)/ai;max=fabs(u0-lastu0);for(i=1;i501;i+)double tran=fabs(ui-lastui);if(maxtran) max=tran;if(max=ep) break;lastbeta=beta;beta=0;for(i=0;i501;i+)beta+=yi*ui;if(fabs(beta-lastbeta)/fabs(beta)=ep) break;return 1/beta;double det(double a501) /求det double det_A=1;doubl

13、e 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;四 程序结果（1）在反幂法中使用三角分解法，且使用的初始迭代向量为u=(1,1,1)时程序运行结果如下所示：(2)在反幂法中使用jacobi迭代法的程序运行结果如下，部分特征值无法

14、计算，部分特征值计算结果正确。五 计算过程中的现象（1）在反幂法程序中，对方程组Auk=yk-1的求解用了三角分解法和JACOBI迭代法两种方法。使用三角分解法时结果正常，但使用JACOBI迭代法时，题干中求解与k（k=1，2，39）最接近的特征值这一问中，部分特征值计算不出来，部分特征值可计算出结果，且结果与使用三角分解法的结果一致，分析原因为部分构造的矩阵不满足JACOBI迭代方法的收敛条件。（2）迭代初始向量的选择对计算结果会产生一定影响，主要表现在收敛速度上。选用初始迭代向量u=(1,0,0,，0)时，程序的运行结果如下所示，可以看出，程序运行结果与使用初始向量u=(1,1,1)的结果差别较大，而且迭代次数明显增加，前者迭代次数最多可到356次，后者最多为129次。使用MATLAB进行计算发现使用u=(1,1,1)向量更加接近准确值。1、初始迭代向量中元素等于0的个数越少，收敛结果越稳定；初始迭代向量中元素等于0的个数越多，收敛结果越不稳定； 2、迭代初始值ui=s(i=1,2,501)且s的绝对值值极大，收敛结果可以稳定但收敛速度减慢，其原因为s的数量级