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

北航数值分析计算实习报告一Companynumber：【0089WT-8898YT-W8CCB-BUUT-202108】北京航空航天大学《数值分析》计算实习报告第一大题学院：自动化科学与电气工程学院专业：控制科学与工程学生姓名：学号：教师：电话：完成日期：2015年11月6日北京航空航天大学BeijingUniversityofAeronauticsandAstronautics实习题目：第一题设有501501的实对称矩阵A，其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0cbieiiaii。矩阵A的特征值为)501,,2,1(ii，并且有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是不同初始向量的情况下的得到的部分结果。（实验结果截图见附录）迭代次数迭代次数159381160381535573350592674611311611304611343611343611343611表1不同初始向量对应的1和501及其迭代次数由表1可以得到如下结论：1.不同的初始向量对本程序的1影响大，对501没有影响，都能保证收敛到正确值。2.初始向量中必须保证501462~uu中至少有一个为1才能保证1收敛到正确值。3.初始向量非零值的多少和大小对迭代次数并没有明显影响。4.为解决初始向量对程序的影响，可以先对A做平移变换再求1。四、实验程序#include#include#includestaticdoubleb=,c=;#definePrecision1e-12voidcopy(doubleb[501],doubley[501]);doubledianji(doublex[],doubley[]);2e,迭代次数为：%d\n,,;printf(最大特征值为：%.12e,迭代次数为：%d\n,,;/**********************************//*以下利用反幂法求解模最小的特征值*//**********************************/InitC(C);2e,迭代次数为：%d\n,,;printf(A的行列式为：%.12e；,;printf(A的条件数cond(A)2为：%.12e\n,;/**********************************************//*以下利用反幂法求与数列Uk中元素最相近的特征值*//**********************************************/doubleUk[39];//保存Uk的值并保存与Uk[i]最接近的λdoubleB[39];intdiedai;iterations=&diedai;for(i=1;i=39;i++){Uk[i-1]=+i*-/40;InitC(C);Initu(u);for(intj=0;j501;j++)C[2][j]-=Uk[i-1];DoolittleC(C,501,2,2);B[i-1]=Get_min_Eigenvalue(C,u,501,2,2,iterations);B[i-1]+=Uk[i-1];printf(μ%-2d=%,与μ%-2d最相近的特征值λ=%\n,i,Uk[i-1],i,B[i-1]);}return0;}doubledetA(doubleC[5][501]){inti=0;doublee=1;for(i=0;i501;i++)e*=C[2][i];returne;}voidInitMatrix(double*p){for(inti=1;i=501;i++){*p=p++;}}voidInitu(double*p){for(inti=1;i=501;i++){if(i=500)*p=0;else*p=1;p++;}//*p=1;}voidA_sub_minI(double*a,doublemin){for(inti=1;i=501;i++){*a=*a-min;a++;}}doubleNeiJi(double*a,double*b){doublee=;for(inti=0;i501;i++){e=e+(*a)*(*b);a++;b++;}returne;}voidget_y(double*y,double*u){doubleNorm_u=sqrt(NeiJi(u,u));for(inti=0;i501;i++){*y=*u/Norm_u;y++;u++;}}voidget_u(double*u,double*y,double*a){u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3];for(inti=2;i499;i++)u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];u[499]=c*y[497]+b*y[498]+a[499]*y[499]+b*y[500];u[500]=c*y[498]+b*y[499]+a[500]*y[500];}voidInitC(doubleC[5][501]){inti;for(i=2;i501;i++){C[0][i]=;C[4][i-2]=;}for(i=1;i501;i++){C[1][i]=;C[3][i-1]=;}for(i=1;i=501;i++){C[2][i-1]=}}doubleGet_Fabs_Eigenvalue(double*a,double*u,int*iterations){doubley[501],B_k0=0,B_k1=0;doublewucha;inti=0;while(1){++i;get_y(y,u);get_u(u,y,a);B_k1=NeiJi(y,u);//是否判断B_K1是否为0wucha=(fabs(B_k1-B_k0))/(fabs(B_k1));//getwuchaif(wucha=Precision)break;elseif(i10000){printf(迭代次数超长，请更改初始向量\n);break;}elseB_k0=B_k1;}*iterations=i;returnB_k1;}doubleGet_min_Eigenvalue(doubleC[5][501],double*u,intn,ints,intr,int*iterations){doubley[501]={0},B_k0=0,B_k1=0;doubleb[501]={0};//储存y值的中间向量doublewucha;inti=0;while(1){++i;g