北航数值分析大作业全.docVIP

学号姓名

学号

姓名

指导老师

《数值分析》计算实习报告

第一大题

题目要求

已知501*501阶的带状矩阵A，其特征值满足...。试求：

，和的值；

A的与数最接近的特征值（k=1，2，...，39）；

A的（谱范数）条件数和行列式detA。

算法设计方案

题目所给的矩阵阶数过大，必须经过去零压缩后进行存储和运算，本算法中压缩后的矩阵如下所示。

由矩阵A的特征值满足的条件可知与之间必有一个最大，则采用幂法求出的一个特征值必为其中的一个：当所求得的特征值为正数，则为；否则为。在求得与其中的一个后，采用带位移的幂法则可求出它们中的另一个，且位移量即为先求出的特征值的值。用反幂法求得的特征值必为。由条件数的性质可得，为模最大的特征值与模最小的特征值之比的模，因此，求出，和的值后，则可以求得。

对于题目的第二问。当我们将位移量设为，再用反幂法求得的相应的特征值即为与之最接近的特征值。

由于在使用反幂法或带位移的反幂法过程中需要解线性方程组，对此我采用三角分解法先分解矩阵，再求迭代所需的解向量。当把进行三角分解后，三角矩阵的对角线元素之积即为detA。

注：算法中涉及求范数时，均为。本报告的非程序部分的矩阵和向量下标均以1开始，程序部分的下标按程序语法规定从0开始计。

创新点

在幂法求按模最大的特征值过程中，归一化后得到再进行迭代求，迭代公式为

为A的非0元素与对应的的积的和，而在程序中变为（为A压缩后的矩阵，其中为501×1阶）

以第一、二、三列的乘积为列，

由上可以看出前三列的相乘在编程时是不太方便的，为了统一和简化编程，做出如下两处改动即可解决问题：

首先，为了方便压缩矩阵的赋值，将的0元素处理后变为如下的矩阵，

；

其次，将向量的空间拓展为505×1，并且有

这样就可以将这个迭代部分统一为

，j=1,2,...,501

四、运算结果及分析

以下计算结果截图为初始迭代向量的情况下算得的。

本算法采用的是各个分量相同且为1的初始迭代向量，运算后得到了正确的结果。将各个分量全改为0后无法运算，但改为相同的其他实数后也能得到相同的正确结论，这是运算时间开销不同。当模最大或最小特征值对应的初始迭代向量的分量为0而其他分量不全为0的时候也能够得到相同的正确结果。原因为即使所求的特征值对应的初始迭代向量分量为0，那么由于计算过程有舍入误差的影响，必然会在迭代的某一步产生一个该分量不为0的新迭代向量。因此，也能得到正确的结论。

五、源程序

#includeiostream.h

#includemath.h

#includeiomanip.h

#defineCM5

#defineCN501

#defineG0.16//矩阵A中的b

#defineH-0.064//矩阵A中的c

voidinitial(doublea[CM][CN]);

intmax(inta,intb);

intmin(inta,intb);

doublepower(double(*a)[CN],doublep);

voidbreakLU(double(*a)[CN]);

voidfun(double(*a)[CN],doubleu[CN],doublex[CN]);

doubledet(double(*a)[CN]);

doublerepower(double(*a)[CN],doublep);

voidmain()

{

doublea[CM][CN],d,t,m,n,s,max,min;

inti;

initial(a);//初始化矩阵A

breakLU(a);

d=det(a);//求行列式值并输出

cout矩阵A的行列式值detA=setprecision(12)setiosflags(ios::scientific)d

endl;

coutendl;

coutendl;

initial(a);//初始化矩阵A

m=power(a,0);//幂法求按模最大的特征值

if(m=0)

{

max=m;

min=power(a,m);//带位移的幂法求另一个特征值

}

else

{

min=m;

max=power(a,m);//带位移的幂法求另一个特征值

}

cout第1个特征值λ1=setprecision(12)setiosflags(ios::scientific)minendl;

cout第501个特征值λ501=setprecision(12)setiosflags(ios::scientific)max

endl;

initial(a);//初始化矩阵A

s=repower(a,0);//反幂法求按模最小的特征值

cout第s个特征值,即