《数值分析》上机实习.docx

1、数值分析上机实习2006级研究生数值分析上机试题 院系 电气工程学院 专业 测检技术与自动化装置 姓名 姜启 学号 2006020067第一题题目：已知12.384120 2.115237 -1.061074 1.112336 -0.113584 0.718719 1.742382 3.067813 -2.0317432.115237 19.141823 -3.125432 -1.012345 2.189736 1.563849 -0.784165 1.112348 3.123124-1.061074 -3.125432 15.567914 3.123848 2.031454 1.836742

2、 -1.056781 0.336993 -1.010103A= 1.112336 -1.012345 3.123848 27.108437 4.101011 -3.741856 2.101023 -0.718280 -0.037585-0.113584 2.189736 2.031454 4.101011 19.897919 0.431637 -3.111223 2.121314 1.7843170.718719 1.563849 1.836742 -3.741856 0.431637 9.789365 -0.103458 -1.103456 0.2384171.742382 -0.78416

3、5 -1.056781 2.101023 -3.111223 -0.103458 14.713846 3.123789 -2.2134743.067813 1.112348 0.336993 -0.718280 2.121314 -1.103456 3.123789 30.719334 4.446782-2.031743 3.123124 -1.010103 -0.037585 1.784317 0.238417 -2.213474 4.446782 40.00001b=(2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343

4、,1.1101230,4.719345,-5.6784392); （1） 用Hoseholder变换，把A化为三对角矩阵B（并打印B）（2） 用超松驰法求解BX=b（取松驰因子W=1.4，X0=0，迭代九次）（3） 用列主元素消去法求解BX=b解题方法的理论依据 ：Household变换：任何实非奇异矩阵A，都可通过初等反射阵相似变化化为HOSSENBERG阵，因此若A是对称阵，则可化为相似的三对角阵，从而解决了实对称三对角化问题。本题中，可依次对矩阵的列向量应用Household变换，完成矩阵的三对角化。算法：（1） 令A0=A, aij(1)=aij,已知Ar-1即Ar-1=(aij(r)

5、；（2） Sr=(ni=r+1(air (r) 2)1/2 ;（3） ar=Sr2+|a(r)r+1,r|Sr;ur=0,0,a(r)r+1,r+Sign(a(r)r+1,r)Sr, a(r)r+2,r, ,a(r)nrT（4） yr=Ar-1ur/ar ;（5） kr=0.5*urTyr/ar ;（6） qr=yr-krur ;Ar=Ar-1-(qrurT+urqrT) ; r=1,2.n-2 .松驰法：在GS方法法已求出x(m)基础上，经过重新组合而得的新序列，此新序列使收敛加快算法如下：列列主消元法：为改善Guass消去法的缺陷主程序如下：#include #include main()

6、int n=9,r,i,j,k; float sr=0,h=0,xr,kr=0; float w=1.4,c=0,c1=0,c2=0,c3=0,c4=0,c5=0; float z3,d9; float ur91; float urt19; float yr91; float qr91; float qrt19; float m99; float t99; float x910=0; float g9=2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392; floa

7、t a99=12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743, 2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124, -1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103, 1.112336,-1.012345,3.123848,27.10843

8、7,4.101011,-3.741856,2.101023,-0.71828,-0.037585, -0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317, 0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417, 1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2

9、.213474, 3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782, -2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001;/*HOUSEHLODER 变换*/for(r=0;r7;r+) /*总循环次数，r=n-2 */sr=0;for(i=r+1;i9;i+) /*求SR*/sr=air*air+sr; sr=sqrt(sr); xr=sr*sr+fabs(ar

10、+1r)*sr; /*求公式中的ar*/for(i=0;i=r;i+) /*求向量ur，把ur设成一列九行的数组*/uri0=0; urr+10=ar+1r+sr*ar+1r/fabs(ar+1r); /*符号函数用A/|A|的形式表式*/for(i=r+2;i9;i+)uri0=air;for(j=0;j9;j+) /*求向量yr*/ h=0; for(i=0;i9;i+) /*求矩阵Ar-1与ur的乘积*/ h=aji*uri0+h; h=h/xr; yrj0=h;for(i=0;i9;i+) /*求kr*/urt0i=uri0;kr=0;for(i=0;i9;i+)kr=urt0i*yr

11、i0+kr; kr=kr/(2*xr);for(i=0;i9;i+) /*求qr*/qri0=yri0-kr*uri0;for(i=0;i9;i+)qrt0i=qri0;for(j=0;j9;j+) /*求qr乘ur的转置*/ for(i=0;i9;i+) mji=qrj0*urt0i;for(j=0;j9;j+) /*求ur乘qr的转置*/ for(i=0;i9;i+) tji=urj0*qrt0i;for(j=0;j9;j+) /*把上面求得的两个矩阵相加*/ for(i=0;i9;i+) tji=mji+tji;for(j=0;j9;j+) /*求Ar*/ for(i=0;i9;i+)

12、aji=aji-tji;printf(houshold is :n);for(j=0;j9;j+) for(i=0;i9;i+) printf(%.6f ,aji); printf( n);/*超松驰法求方解*/ for(i=1;i10;i+) /*迭代九次*/ for(j=1;j9;j+) c2=c2+a0i*xji-1/a00; /*求X0的一次迭代*/ c=w*(g0/a00-c2); x0i=c-0.4*x0i-1; c2=0; for(j=1;j8;j+) /*求X1到X7的一次迭代*/ for(k=0;kj;k+) c1=c1+ajk*xki/ajj; /*公式中的第一个求和式*/

13、 for(k=j+1;k9;k+) c2=c2+ajk*xki-1/ajj; /*公式中的第二个求和式*/ c=w*(gj/ajj-c1-c2); xji=c-0.4*xji-1; c1=0; c2=0; for(j=0;j8;j+) c1=c1+a8j*xji/a88; /*求X8的一次迭代*/ c=w*(g8/a88-c1); x8i=c-0.4*x8i-1; c1=0; printf(sor result are:n); for(i=0;i9;i+) printf(%.5f, ,xi8); printf(n);/*列主元素求解*/ for(i=0;iaii) for(j=i,k=0;ji+3;j+,k+) zk=aij;aij=ai+1j;ai+1j=zk; c3=gi;gi=gi+1;gi+1=c3; c5=ai+1i; for(j=i;ji+3;j+) ai+1j=ai+1j-aij*c5/aii; gi+1=g