1、北京航空航天大學數值分析大作業一學院名稱自動化專業方向控制工程學 號 ZY1403140學生姓名許陽教 師 孫玉泉日 期2021年11月26日設有501 501的實對稱矩陣A,a1b c bA ccbc ba50i0.1其中，ai (1.64 0.024i)sin(0.2i) 0.64e-r(i 1,2, ,501),b 0.16,c0.064。矩陣A的特征值為i(i 1,2, ,501),并且有12501,| s| 2%| i|1,501和s的值。A的與數k 1 k10-1最接近的特征值i (k 1,2, ,39)。40kA的(譜范數)條件數8口(人)2和行列式detAo方案設計1求1 ,5

2、01和s的值。s為按模最小特征值，| s| min | i |0可使用反幕法求得。1 i 5011, 501分別為最大特征值及最小特征值。可使用幕法求出按模最大特征值，如結果為正，即為501,結果為負，那么為10使用位移的方式求得另一特 征值即可。2求A的與數k 1 k二0-1最接近的特征值i (k 1,2,.,39) o 40k題目可看成求以k為偏移量后，按模最小的特征值。即以 k為偏移量做位移，使用反幕法求出按模最小特征值后，加上即為所求。3求A的(譜范數)條件數8門/)2和行列式detAo矩陣A為非奇異對稱矩陣，可知，(1-1)cond(A)2 |max | min頁腳下載后可刪除，如有

3、侵權請告知刪除!min 為按模最小特征值。其中 max 為按模最大特征值，(2-1)(2-2) 那么(2-3)(2-4)detA 可由 LU 分解得到。因 LU 均為三角陣，那么其主對角線乘積即為 A 的行列式。二 算法實現1 冪法使用如下迭代格式：任取非零向量u0(u1(0), , un(0)Tyk 1 uk 1/max |uk 1 | uk Ayk 1k sgn(max |uk 1 |) max|uk 1 |終止迭代的控制理論使用 | k k 1 |/ | k |，實際使用| k | | k1|/| k |由于不保存A矩陣中的零元素，只保存主對角元素a501及b,c值上式中 uk Ayk

4、1簡化為：u(1) a(1)y(1) by(2) cy(3)u(2) by(1) a(2)y(2) by(3) cy(4)u(500) cy(498) by(499) a(500)y(500) by(501) u(501) cy(499) by(500) a(501)y(501) u(i) cy(i 2) by(i 1) a(i)y(i) by(i 1) cy(i 2) (i 3, ,499)2 反冪法使用如下迭代格式：任取非零向量u0(h1(0), ,hn(0)Tyk 1 uk 1/max|uk 1 | -1 uk A yk 1ksgn(max | uk 1 |) max |uk 1 |其中

5、 uk A 1yk 1Aukyk 1 ，解方程求出 uk 。頁腳下載后可刪除，如有侵權請告知刪除！求解過程中使用LU分解，由于A為5對角矩陣，選擇追趕法求取LU分解求解過程如下：LUuk yk 1Lxk yk 1Uuk xk uk頁腳下載后可刪除，如有侵權請告知刪除!追趕法求LU分解的實現:a1 b cbA cc LUb(2-5)c ba501p11 t1 q1r2z3q499t 500z501r501p5011由上式推出分解公式如下：p1a1,t1b/a1r2b, p2a2 r2t1qic/ pi,i1,.,499ti (b riqi 1)/ pi,i 2,.,500(2-6)zic,i 3

6、,.,501rib cti 2,i 3,.,501piai cqi 2 riti 1,i 3,.,501推導出回代求解公式如下：x1y1 / p1x2(y2 r2x1) / p2(2-7)xi(yizixi 2rixi 1)/ pi ,i3,.,501U501 X50iU500 X500 t500u501UiXi tiUi 1 qiXi 2,i 499，,13 cond(A)2及A行列式求解| i |cond (A)2 1s |由式(2-5)可得：501 det A pi i 1(2-8)(2-9)(2-10)三源程序#include <stdio.h>#include <m

7、ath.h>double ep=1e-12,b=0.16,c=-0.064;int j=0;double power(double a501); 幕法double inv_power(double a501); /反幕法double det(double a501) ;/求 detint main()主程序int i,k;double A501,B501,beta_1,beta_501,beta_s,beta_k;double mu;for(i=0;i<501;i+)Ai=(1.64-0.024*(i+1)*sin(0.2*(i+1)-0.64*exp(0.1/(i+1);beta

8、_1=power(A);第一問printf(" t= %.12et 迭代次數：%dn",beta_1,j);for(i=0;i<501;i+)/位移Bi=Ai-beta_1;beta_501=power(B)+beta_1;printf(" 入 501 %.12et 一代次數:%dn",beta_501,j);beta_s=inv_power(A);printf(" t= %.12et 迭代次數:%dn",beta_s,j);for(k=1;k<=39;k+)第二問mu=beta_1+k*(beta_501-beta_1)

9、/40;for(i=0;i<501;i+)Bi=Ai-mu;beta_k=inv_power(B)+mu;printf("入 i%.12et 迭代次數：dn”,k,beta_k,j);printf("cond(A)2= %.12en",beta_1/beta_s);/第三問printf("detAt= %.12en",det(A);double power(double a501)/冪法int i=0,N=5000;double b=0.16,c=-0.064;double u501,y501;double m=1,beta;for(i=

10、0;i<501;i+)ui=1;j=0;while(j<N)for(i=0;i<501;i+)yi=ui/fabs(m);u0=a0*y0+b*y1+c*y2;u1=b*y0+a1*y1+b*y2+c*y3;u499=c*y497+b*y498+a499*y499+b*y500;u500=c*y498+b*y499+a500*y500;for(i=2;i<499;i+)ui=c*yi-2+b*yi-1+ai*yi+b*yi+1+c*yi+2;beta=0;for(i=0;i<501;i+)if(fabs(ui)>=fabs(beta) beta=ui;if(

11、beta<0)if(fabs(fabs(beta)-fabs(m)/fabs(beta)<ep) break;if(fabs(beta-m)/fabs(beta)<ep) break;m=beta;j+;return beta;double inv_power(double a501)/反冪法double p501,r501,t501,q501,u501,y501;double beta,m=1;int i,N=1000;p0=a0;t0=b/p0;r1=b;p1=a1-r1*t0;q0=c/p0;q1=c/p1;t1=(b-r1*q0)/p1;for(i=2;i<50

12、1;i+)ri=b-c*ti-2;pi=ai-c*qi-2-ri*ti-1;qi=c/pi;ti=(b-ri*qi-1)/pi;for(i=0;i<501;i+)ui=1;j=0;while(j<N)for(i=0;i<501;i+)yi=ui/fabs(m);u0=y0/p0;u1=(y1-r1*u0)/p1;for(i=2;i<501;i+)ui=(yi-c*ui-2-ri*ui-1)/pi;u499=u499-t499*u500;for(i=498;i>=0;i-)ui=ui-ti*ui+1-qi*ui+2;beta=0;for(i=0;i<501;i

13、+)if(fabs(ui)>=fabs(beta)beta=ui;if(beta<0)if(fabs(fabs(beta)-fabs(m)/fabs(beta)<ep) break;if(fabs(beta-m)/fabs(beta)<ep)break;m=beta;j+;return 1/beta;double det(double a501) / 求 detdouble det_A=1;double p501,r501,t501,q501;int i;p0=a0;t0=b/p0;r1=b;p1=a1-r1*t0;q0=c/p0;q1=c/p1;t1=(b-r1*q0

14、)/p1;for(i=2;i<501;i+)ri=b-c*ti-2;pi=ai-c*qi-2-ri*ti-1;qi=c/pi;ti=(b-ri*qi-1)/pi;for(i=0;i<501;i+)det_A=det_A*pi;return det_A;頁腳下載后可刪除，如有侵權請告知刪除！四程序結果'&Progrsm Files ：>8& C-Fr« b CJvJfc:Kfinall.exf"人工=-1.07Q011361488e*061入501=9.72463 41 Bl479e*00Xs=-5.S5791W794214R-HB3

15、X ii-1 .01629340331601XiZ-?.5#5707425B68C*000Xi3=-9.172&7W3V28e*tiUl!lX i4-8.65228007898e*000X 15=-8.09348380O«Xit-?*G 5?4日5407G摩甘也。目Xi?=-7 *工1 yBJMtHHbTL c *修婭>Xi8-£ .6il76O39397»««eaAi9-6 _抽帚 以2瞿6耳9氫4mMmXue- 51。工的2G28。*0O0Am-5 .lld083529#12c 現加Xii2-4-578872176865e0O

16、0Xii3=-4.07»2«fi7*MA +HflB1 ild-3 ,ESd2112157Elo *G(J0Xii5-3國41町加上8133H即90A lib=T,b2b4W4byi3«e*MUMA117=-2.0032307t95G4e，的01 119=-l-E03&E7fcli227o«QO0Xii9-9.935&O606Sa?Se USiki2b=-4 _87B42f>73B85We -HBlA 121=2.2317362495?5e-002X 152=J 蒐妙一曲 1襖次-W次次次.次次次次.次次次次次砍次次次次次次次以 弋

17、弋弋弋弋弋弋弋弋弋弋弋弋弋弋弋弋弋弋七弋4>心弋弋 中勺號勺寸zfEKr= ff m1 1母 1S 1 "I，71. I F? x-TAh l ,- Ta Ta IIT1 111111.1 TA - A t-ta - Ta -Bl -.i -Jr- .1、工11178fi2226 slfi-lt-7b2712 4 14- 644- 0Rdetft = 2.772'4f b1417b2e +118清報任意褥緡續,.一 _uuniKA；2- 1.92520427;B83e*003*tfcProgram Files xfiGAC-free 5 CJVJ®|Ri

18、63;finall.eKe'回345&7UO90 1234 5 67R9 2222222233323333 2L i £ 工 L L 七工 工 t -L al 七 i .n*1 工 xx入AxkxxxxA 篇xkx源117b 4 w 4 2 2 9 4 ft 413 121 次慶次次次次:A次火-A/.外次次次攻收 七七弋弋tt七七七弋弋弋t弋弋 . TA-1- - Ta - a rl* 1 i 工、工 .I、工.工 11.1、Ta%-Ta-1 4 052 8?8?626?3e +000 -1.589445881881。*閱 -2演網8水口273*鬧口 -2 4 5E#0r?£S97073* *000 -3.08O2BS0S307e *600 -3. bH«7?B'/aS4SMe +0WM -4.0913?85104Sle+00G =4.37S2?9e-5.13292