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

整理文档很辛苦,赏杯茶钱您下走!

免费阅读已结束,点击下载阅读编辑剩下 ...

阅读已结束,您可以下载文档离线阅读编辑

资源描述

北京航空航天大学《数值分析》计算实习报告第一大题学院:自动化科学与电气工程学院专业:控制科学与工程学生姓名:学号:教师:电话:完成日期: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

1 / 26
下载文档,编辑使用

©2015-2020 m.777doc.com 三七文档.

备案号:鲁ICP备2024069028号-1 客服联系 QQ:2149211541

×
保存成功