北航大学数值分析计算实习报告.doc

北航大学数值分析计算实习报告 实习题目： 第一题 设有的实对称矩阵A， 其中，。矩阵A的特征值为，并且有 1.求，和的值。 2.求A的与数最接近的特征值。 3.求A的(谱范数)条件数和行列式detA。 说明： 在所用的算法中，凡是要给出精度水平的，都取。 选择算法时，应使矩阵A的所有零元素都不储存。 打印以下内容： 全部源程序； 特征值以及的值。 4.采用e型输出实型数，并且至少显示12位有效数字。 一、算法设计方案 1、求，和的值。 由于，可知绝对值最大特征值必为和其中之一，故可用幂法求出绝对值最大的特征值，如果=0，则=，否则=。将矩阵A进行一下平移： （1） 对用幂法求出其绝对值最大的特征值，则A的另一端点特征值或为+。 为按模最小特征值，，可对A使用反幂法求得。 2、求A的与数最接近的特征值。 计算-，其模值最小的值对应的特征值与最接近。因此对A进行平移变换： （2） 对用反幂法求得其模最小的特征值，则=+。 3、求A的(谱范数)条件数和行列式detA。 由矩阵A为非奇异对称矩阵可得： （3） 其中为按模最大特征值，为按模最小特征值，通过第一问我们求得的和可以很容易求得A的条件数。 在进行反幂法求解时，要对A进行LU分解得到。因L为单位下三角阵，行列式为1，U为上三角阵，行列式为主对角线乘积，所以A的行列式等于U的行列式，为U的主对角线的乘积。 二、 算法实现 矩阵存储 原矩阵A为一个上、下半带宽都为2的501×501的带状矩阵，由于矩阵中的0元素太多，如果分配一个501×501的空间保存矩阵的话会浪费很多空间。因此，为了节省存储量，A的带外元素不给存储，值存储带内元素，如下C矩阵所示： （4） C是一个5×501的矩阵，相比A大大节省了存储空间，在数组C中检索矩阵A的带内元素的方法是： （5） 幂法 幂法迭代公式如下： （6） 其中，不断迭代当时即可认为其满足精度要求，令。 在程序中计算时，根据A矩阵的特点，简化如下： （7） 3、反幂法 反幂法迭代公式如下： （8） 当k足够大时，。 在求解时，可先对A进行Doolittle分解，由于A是带状结构，所以分解出的L、U也是带状结构，利用C矩阵进行Doolittle分解并求向量的算法如下： （1）作分解A=LU 对于执行： （9） 由于C语言中数组下标是从0开始的，所以在程序中矩阵元素c的下标都减1。 求解（数组b先是存放原方程组右端向量，后来存放中间向量y,在程序中b和y都保存在数组y[501]中。） （10） 求出后，其他部分与幂法求解相同。  结果分析 实验表明，本程序中，初始向量对结果影响较大，合适的初始向量对得到正确的收敛结果比较重要，如表1是不同初始向量的情况下的得到的部分结果。（实验结果截图见附录） 迭代次数 迭代次数 159 381 160 381 535 573 350 592 674 611 311 611 304 611 343 611 343 611 343 611 表1 不同初始向量对应的和及其迭代次数 由表1可以得到如下结论： 不同的初始向量对本程序的影响大，对没有影响，都能保证收敛到正确值。 初始向量中必须保证中至少有一个为1才能保证收敛到正确值。 初始向量非零值的多少和大小对迭代次数并没有明显影响。 为解决初始向量对程序的影响，可以先对A做平移变换再求。 实验程序 #include stdio.h #include math.h #include stdlib.h static double b=0.16,c=-0.064; #define Precision 1e-12 void copy(double b[501],double y[501]); double dianji(double x[],double y[]);//计算两个向量内积 void InitMatrix(double *p);//Init Matrix A double NeiJi(double a[],double b[]);//get 2范数 void get_y(double *y,double *u);//get y void get_u(double *u,doubl