《北航数值分析计算实习报告一()》由会员分享，可在线阅读，更多相关《北航数值分析计算实习报告一()（22页珍藏版）》请在金锄头文库上搜索。
1、北京航空航天大学数值分析计算实习报告第一大题学 院：自动化科学与电气工程学院专 业： 控制科学与工程 学 生 姓 名： 学 号： 教 师： 电 话： 完 成 日 期： 2015年11月6日 北京航空航天大学Beijing University of Aeronautics and Astronautics实习题目：第一题 设有的实对称矩阵A，其中，。矩阵A的特征值为，并且有1.求，和的值。2.求A的与数最接近的特征值。3.求A的(谱范数)条件数和行列式detA。说明：1. 在所用的算法中，凡是要给出精度水平的，都取。2. 选择算法时，应使矩阵A的所有零元素都不储存。3. 打印以下内容：（1）
2、全部源程序；（2） 特征值以及的值。4.采用e型输出实型数，并且至少显示12位有效数字。一、算法设计方案1、求，和的值。由于，可知绝对值最大特征值必为和其中之一，故可用幂法求出绝对值最大的特征值，如果=0，则=，否则=。将矩阵A进行一下平移： （1）对用幂法求出其绝对值最大的特征值，则A的另一端点特征值或为+。为按模最小特征值，可对A使用反幂法求得。2、求A的与数最接近的特征值。计算-，其模值最小的值对应的特征值与最接近。因此对A进行平移变换： （2）对用反幂法求得其模最小的特征值，则=+。3、求A的(谱范数)条件数和行列式detA。由矩阵A为非奇异对称矩阵可得： （3）其中为按模最大特征值，
3、为按模最小特征值，通过第一问我们求得的和可以很容易求得A的条件数。在进行反幂法求解时，要对A进行LU分解得到。因L为单位下三角阵，行列式为1，U为上三角阵，行列式为主对角线乘积，所以A的行列式等于U的行列式，为U的主对角线的乘积。二、 算法实现1、 矩阵存储原矩阵A为一个上、下半带宽都为2的501501的带状矩阵，由于矩阵中的0元素太多，如果分配一个501501的空间保存矩阵的话会浪费很多空间。因此，为了节省存储量，A的带外元素不给存储，值存储带内元素，如下C矩阵所示： （4）C是一个5501的矩阵，相比A大大节省了存储空间，在数组C中检索矩阵A的带内元素的方法是： （5）2、 幂法 幂法迭代
4、公式如下： （6）其中，不断迭代当时即可认为其满足精度要求，令。在程序中计算时，根据A矩阵的特点，简化如下： （7）3、反幂法反幂法迭代公式如下： （8）当k足够大时，。在求解时，可先对A进行Doolittle分解，由于A是带状结构，所以分解出的L、U也是带状结构，利用C矩阵进行Doolittle分解并求向量的算法如下：（1）作分解A=LU 对于执行： （9） 由于C语言中数组下标是从0开始的，所以在程序中矩阵元素c的下标都减1。（2） 求解（数组b先是存放原方程组右端向量，后来存放中间向量y,在程序中b和y都保存在数组y501中。） （10）求出后，其他部分与幂法求解相同。三、 结果分析实验
5、表明，本程序中，初始向量对结果影响较大，合适的初始向量对得到正确的收敛结果比较重要，如表1是不同初始向量的情况下的得到的部分结果。（实验结果截图见附录）迭代次数迭代次数159381160381535573350592674611311611304611343611343611343611表1 不同初始向量对应的和及其迭代次数由表1可以得到如下结论：1. 不同的初始向量对本程序的影响大，对没有影响，都能保证收敛到正确值。2. 初始向量中必须保证中至少有一个为1才能保证收敛到正确值。3. 初始向量非零值的多少和大小对迭代次数并没有明显影响。4. 为解决初始向量对程序的影响，可以先对A做平移变换再求
6、。四、 实验程序#include #include #include static double b=0.16,c=-0.064;#define Precision 1e-12void copy(double b501,double y501);double dianji(double x,double y);/计算两个向量内积void InitMatrix(double *p);/Init Matrix Adouble NeiJi(double a,double b);/get 2范数void get_y(double *y,double *u);/get yvoid get_u(doubl
7、e *u,double *y,double *a);/get uvoid Initu(double *p);/初始化初始向量udouble Get_Fabs_Eigenvalue(double *a,double *u,int *iterations);/循环迭代得到绝对值最大特征值void A_sub_minI(double *a,double min);/A-min*Ivoid InitC(double C5501);/初始化数组Cvoid DoolittleC(double C5501,int n,int s,int r);/进行Doolittle分解int min(int a,int
8、b);/取返回a,b最小值int max(int a,int b);/求最大值void Doolittle_getx(double C5501,double y501,double u501,int n,int s,int r);double Get_min_Eigenvalue(double C5501,double *u,int n,int s,int r,int *iterations);double detA(double C5501);/求A的行列式struct FinalValuedouble min;/特征值最小值double max;/特征值最大值double abs_min;
9、/模最小特征值double abs_max;/模最大特征值double detA;/A的行列式double cond2;/A的条件数int min_iterations;/最小值迭代次数int max_iterations;/最大值迭代次数int abs_min_iterations;/求绝对值最小的迭代次数;int main()FinalValue main_num= 0,0,0,0;double temp;/两值交换中间变量int temp1;double u501=0;/,y501;/为u0赋初值double Norm_u=0;/范数double C5501 = 0;InitC(C);/
10、初始化Cint i=0,*iterations;Initu(u);iterations = &main_num.min_iterations;main_num.min = Get_Fabs_Eigenvalue(C2,u,iterations);/将求得的绝对值最大特征值放到min变量中main_num.abs_max = main_num.min;A_sub_minI(C2,main_num.min);for(i=0;imain_num.max)temp = main_num.min;main_num.min = main_num.max;main_num.max = temp;temp1
11、= main_num.min_iterations;main_num.min_iterations = main_num.max_iterations;main_num.max_iterations = temp1;printf(最小特征值为：%.12e,迭代次数为：%dn,main_num.min,main_num.min_iterations);printf(最大特征值为：%.12e,迭代次数为：%dn,main_num.max,main_num.max_iterations);/*/*以下利用反幂法求解模最小的特征值*/*/InitC(C);/初始化CInitu(u);Doolittle
12、C(C,501,2,2);main_num.detA = detA(C);iterations = &main_num.abs_min_iterations;main_num.abs_min = Get_min_Eigenvalue(C,u,501,2,2,iterations);main_num.cond2 = fabs(main_num.abs_max/main_num.abs_min);printf(绝对值最小特征值为：%.12e,迭代次数为：%dn,main_num.abs_min,main_num.abs_min_iterations);printf(A的行列式为：%.12e；,main_num.detA);printf(A的条件数cond(A)2为：%.12en,main_num.cond2);/*/*以下利用反幂法求与数列Uk中元素最相近的特征值*/*/double Uk39;/保存Uk的值并保存与Uki最
