北航数值分析计算实习报告一

北京航空航天大学《数值分析》计算实习报告第一大题学院：自动化科学与电气工程学院专业：控制科学与工程学生姓名：学号：教师：电话：完成日期：2015年11月6日北京航空航天大学BeijingUniversityofAeronauticsandAstronautics实习题目：第一题设有501501的实对称矩阵A，5011Aabcbccbcba其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0cbieiiaii。矩阵A的特征值为)501,,2,1(ii，并且有||min||,501150121iis1.求1，501和s的值。2.求A的与数4015011kk最接近的特征值)39,,2,1(kki。3.求A的(谱范数)条件数2)A(cond和行列式detA。说明：1.在所用的算法中，凡是要给出精度水平的，都取12-10。2.选择算法时，应使矩阵A的所有零元素都不储存。3.打印以下内容：（1）全部源程序；（2）特征值），，39,...,2,1(,s5011kki以及Adet,)A(cond2的值。4.采用e型输出实型数，并且至少显示12位有效数字。一、算法设计方案1、求1，501和s的值。由于||min||,501150121iis，可知绝对值最大特征值必为1和501其中之一，故可用幂法求出绝对值最大的特征值，如果=0，则1=，否则501=。将矩阵A进行一下平移：I-AA'（1）对'A用幂法求出其绝对值最大的特征值'，则A的另一端点特征值1或501为'+。s为按模最小特征值，||min||5011iis，可对A使用反幂法求得。2、求A的与数4015011kk最接近的特征值)39,...,2,1(kki。计算1)1,2,...,50=(ii-k，其模值最小的值对应的特征值k与k最接近。因此对A进行平移变换：)39,,2,1k-AAkk（I（2）对kA用反幂法求得其模最小的特征值'k，则k='k+k。3、求A的(谱范数)条件数2)(Acond和行列式detA。由矩阵A为非奇异对称矩阵可得：||)(minmax2Acond（3）其中max为按模最大特征值，min为按模最小特征值，通过第一问我们求得的和s可以很容易求得A的条件数。在进行反幂法求解时，要对A进行LU分解得到。因L为单位下三角阵，行列式为1，U为上三角阵，行列式为主对角线乘积，所以A的行列式等于U的行列式，为U的主对角线的乘积。二、算法实现1、矩阵存储原矩阵A为一个上、下半带宽都为2的501×501的带状矩阵，由于矩阵中的0元素太多，如果分配一个501×501的空间保存矩阵的话会浪费很多空间。因此，为了节省存储量，A的带外元素不给存储，值存储带内元素，如下C矩阵所示：000000501500499321ccccbbbbbaaaaaabbbbbccccC（4）C是一个5×501的矩阵，相比A大大节省了存储空间，在数组C中检索矩阵A的带内元素ija的方法是：j1,sj-icaA中的元素的带内元素Cij（5）2、幂法幂法迭代公式如下：kTkkkkkkkuyAyuuyuu111k11211kn0/R任取非零向量（6）其中kklim，不断迭代当||/||1kkk时即可认为其满足精度要求，令k。在程序中计算1ukkAy时，根据A矩阵的特点，简化如下：)499,,3()2()1()()()1()2()()501()500()499()501()501()500()499()498()500()4()3()2()2(C)1()2()3()2()1(C)1(]][3[]501][3[]500][3[]2][3[]1][3[iicyibyiyiCibyicyiuyCbycyubyyCbycyucybyybyucybyyui（7）3、反幂法反幂法迭代公式如下：kTkkkkkkkuyyAuuyuu111k11211kn0/R任取非零向量（8）当k足够大时，k1s。在求解1kkyAu时，可先对A进行Doolittle分解，由于A是带状结构，所以分解出的L、U也是带状结构，利用C矩阵进行Doolittle分解并求向量ku的算法如下：（1）作分解A=LU对于n,,2,1k执行：)],min(,,2,1[/)(:)],min(,,1,[:,11),,1max(,1,1,1,11),,1max(,1,1,1,1nrkkkicccccnskkkjcccckskskritkskttstikskikskiksjrktjsjttstkjsjkjsjk（9）由于C语言中数组下标是从0开始的，所以在程序中矩阵元素c的下标都减1。（2）求解yUxbLy,（数组b先是存放原方程组右端向量，后来存放中间向量y,在程序中b和y都保存在数组y[501]中。）)1,,2,1(/)(/),,3,2(,1),min(1,1,11),1max(,1nnicxcbxcbxnibcbbistnsiittstiiinsnniritttstiii（10）求出ku后，其他部分与幂法求解相同。三、结果分析实验表明，本程序中，初始向量Tuuuu501211501,,,对结果影响较大，合适的初始向量对得到正确的收敛结果比较重要，如表1是不同初始向量的情况下的得到的部分结果。（实验结果截图见附录）1501u1迭代次数501迭代次数11u000370809810853.2-e159000787246340993.9e381121uu000e380809810853.2-160000787246340993.9e38114001uu000962085531594.9-e535000087246340989.9e57314611uu000659539789789.9-e350000407246340988.9e59214621uu001509.0700113611-e674000727246340987.9e6111462u001502.0700113611-e311000727246340987.9e6111481u001501.0700113611-e304000727246340987.9e6111501462uu001502.0700113611-e343000727246340987.9e61115011uu001502.0700113611-e343000727246340987.9e611T501,,2,1001502.0700113611-e343000727246340987.9e611表1不同初始向量对应的1和501及其迭代次数由表1可以得到如下结论：1.不同的初始向量对本程序的1影响大，对501没有影响，都能保证收敛到正确值。2.初始向量中必须保证501462~uu中至少有一个为1才能保证1收敛到正确值。3.初始向量非零值的多少和大小对迭代次数并没有明显影响。4.为解决初始向量对程序的影响，可以先对A做平移变换再求1。四、实验程序#includestdio.h#includemath.h#includestdlib.hstaticdoubleb=0.16,c=-0.064;#definePrecision1e-12voidcopy(doubleb[501],doubley[501]);doubledianji(doublex[],doubley[]);//计算两个向量内积voidInitMatrix(double*p);//InitMatrixAdoubleNeiJi(doublea[],doubleb[]);//get2范数voidget_y(double*y,double*u);//getyvoidget_u(double*u,double*y,double*a);//getuvoidInitu(double*p);//初始化初始向量udoubleGet_Fabs_Eigenvalue(double*a,double*u,int*iterations);//循环迭代得到绝对值最大特征值voidA_sub_minI(double*a,doublemin);//A-min*IvoidInitC(doubleC[5][501]);//初始化数组CvoidDoolittleC(doubleC[5][501],intn,ints,intr);//进行Doolittle分解intmin(inta,intb);//取返回a,b最小值intmax(inta,intb);//求最大值voidDoolittle_getx(doubleC[5][501],doubley[501],doubleu[501],intn,ints,intr);doubleGet_min_Eigenvalue(doubleC[5][501],double*u,intn,ints,intr,int*iterations);doubledetA(doubleC[5][501]);//求A的行列式structFinalValue{doublemin;//特征值最小值doublemax;//特征值最大值doubleabs_min;//模最小特征值doubleabs_max;//模最大特征值doubledetA;//A的行列式doublecond2;//A的条件数intmin_iterations;//最小值迭代次数intmax_iterations;//最大值迭代次数intabs_min_iterations;//求绝对值最小的迭代次数};intmain(){FinalValuemain_num={0,0,0,0};doubletemp;//两值交换中间变量inttemp1;doubleu[501]={0};//,y[501];//为u0赋初值doubleNorm_u=0;//范数doubleC[5][501]={0};InitC(C);//初始化Cinti=0,*iterations;Initu(u);iterations=&main_num.min_iterations;main_num.min=Get_Fabs_Eigenvalue(C[2],u,iterations);//将求得的绝对值最大特征值放到min变量中main_num.abs_max=main_num.min;A_sub_minI(C[2],main_num.min);for(i=0;i501;i++)u[i]=i+1;iterations=&main_num.max_iterations;main_num.max=Get_Fabs_Eigenvalue(C[2],u,iterations);//将求得的绝对值最大特征值放到min变量中main_num.max+=main_num.min;if(main_num.minmain_num.max){temp=main_num.min;main_num.min=main_num.max;main_num.max=temp;temp1=main_num.min_iterations