北航数值分析计算实习二.pdf

1、数值 分析 计算 实习 （二 ） 姓名：张时毓 学号：SY1403531 一、方案设计 1. 用带双步位移的 QR 分 解法 求 A 的全体特征根 先对 A 进行拟上三角化，得到 ( 1) ，然后对 ( 1) 用带双步位移的 QR 分解法进行 QR 分解，逐步得到 A 的全部特征根。算法流程图如图 1 所 示。 2. 求实特征根对应的特征向量 设拟上三角化过程为 ( 1) = 1 其中 为正交阵。 QR 分解迭代 过程为 1 = ( 1) +1 = 1 = 1 2 1 1 1 1 1 2 = ( 1 2 ) 1 ( 1) ( 1 2 ) 记 = 1 2 ，则有 +1 = 1 ( 1) = 1

2、1 = ( ) 1 ( ) 其中 为正交阵。 记Q = ，则有 +1 = 1 其中 Q 为正 交阵。 因此 +1 与 A 相似，设 A 和 +1 的特征值分别为 1 ， 2 和 1 ， 2 ，对应的特 征值分别为x 1 ，x 2 x 和x 1 ，x 2 x ，则有 = ( = 1 ，2 ) x = Qx ( = 1 ，2 ) 实特征根对应的特征向量即为 Q 中 对应的列向量。 3. QR 分解 从 r=1 开始 ，记 ( 1) = ，构造相似变换矩阵 H r ，使 ( +1) = ( )r=n-1 时，有 ( ) = 1 ( 1)其中 R = ( )Q = 1 2 1A = QR 开始 初始化

3、 给定精度er0， 最大迭代次数L 对A进行拟上三角化 得到A (n-1) k=1 m=n A 1 =A (n-1) |a m, m-1(k) | er0 ? |a m-1, m-2(k) | er0 ? k=L ? 带双步位移的 QR分解 k=k+1 得到一个 实特征根 m=m-1 得到一对 共轭特征根 m=m-2 M=2 ? 得到最后一个或两个 特征根 得到全部特征根 成功 超过最大迭代次数 失败 结束 Y Y Y N N N图 1 带双步位移的 QR 分解法 求特征根算法流程图 H r 的构造方法如下： 若 = 0 ，( + 2) ，则取 = ； 否则，取 = (0 0, ( ) ( )

4、 ) e = = + ( ) 0 + ( ) 0 = = (0 0, ( ) ( ) ) = 2 2 2 = 1 0 0 1 若每步计算先算 H r ，再 由 ( +1) = ( )计算 ( +1) ， 则需做一次矩阵减法、 一次向量 乘法和一次矩阵乘法，即 n 2 次减法和 n 2 +n 3 次乘法，计算量大，因此，令 = 1 2 2 = = ( ) = ( ) 则有 ( +1) = ( ) ( +1) = ( ) 此时，计算 ( +1) 只需作 n 2 次减法和 n+ n 2 次乘法，计算量大大降低。 4. 拟上三角化 用 QR 方法求 特征值时， 直 接对 A 进行分解迭代次数太多， 因

5、此 先对 A 进行拟上三角化， 再对得到的拟上三角阵进行 QR 分解 ，从而减少迭代次数。 拟上三角化过程与 QR 分 解类似，从 r=1 开始，记 ( 1) = ，构造相似变换矩阵 H r ，使 ( +1) = ( ) r=n-2 时，有 ( 1) = 2 ( 2) 2得到的 ( 1) 即为拟上三角阵。 H r 的构造方法如下： 若 = 0 ，( + 1) ，则取 = ； 否则，取 = (0 0, +1, ( ) ( ) ) e = +1 = + +1, ( ) 0 + +1, ( ) 0 = +1 = (0 0, ( ) ( ) ) = 2 2 2 = 0 0 为减小计算量，令 = 1 2

6、 2 = = ( ) = ( ) = = 则有 ( +1) = ( ) 5. 基本 QR 方法 令 1 = ，对 A 作 QR 分解 = +1 = = 1 +1 相似于 ， 即特征值相同。 产生的矩阵序列 收敛于分块上三角阵， 其对角块均 为一阶和二阶子块， 且对角块中每个一阶子块给出 A 的一个实特征值， 每个二阶子块给出 A 的一对复共轭特征值。 6. 带双步位移的 QR 方法 为加快收敛速度，对基本 QR 分解法 进行修正，提出带双步位移的 QR 方法 。 令 = 2 1 ( ) + 2 ( ) + 1 ( ) 2 ( ) 对 进行 QR 分解 = 则有 +1 = 与基本 QR 法相比，

7、每次迭代过程虽然增加了两次矩阵乘法，但减少一次 QR 分 解，且 避免了复数运算。 二、源代码 #include #include #define dim 10 /矩阵维数 #define er0 1E-12 /精度 #define L 1000 /最大迭代次 数 double Adimdim=0,0; double Indimdim=0,0; /单位阵 double FFdimdim=0,0; /缓存数据 double Qdimdim=0,0; /相似变换矩 阵 A 分块对角阵 double Qfdimdim=0,0; void display_u(double xdim) / 向量输出函

8、数 long i=0; long j=0; for(i=0;idim;i+) printf(“%12.12et“,xi); printf(“n“); void display_m(double aadimdim) / 矩阵 输出函数 long i=0; long j=0; for(i=0;idim;i+) for(j=0;jdim;j+) printf(“%12.12et“,aaij); printf(“n“); printf(“n“); double length(double xdim) / 计算向量的模 double leng=0; long i=0; for(i=0;idim;i+)

9、leng=leng+xi*xi; leng=sqrt(leng); return leng; double vtr_vtr(double x1dim,double x2dim) /1*n 向 量与 n*1 向量相乘 double vv=0; long i=0; for(i=0;idim;i+) vv=vv+x1i*x2i; return vv; void mtx_vtr(double aadimdim,double xdim,double mvdim) /n*n 矩阵与 n*1 向量相乘 /double mvdim=0; long i=0,j=0; for(i=0;idim;i+) mvi=0

10、; for(j=0;jdim;j+) mvi=mvi+aaij*xj; void vtr_vtr2(double x1dim,double x2dim,double mmdimdim) /n*1 向 量与 1*n 向量 相乘 long i=0,j=0; for(i=0;idim;i+) for(j=0;jdim;j+) mmij=x1i*x2j; void vtr_sub(double x1dim,double x2dim,double uudim) /向量- 向量 long i=0; for(i=0;idim;i+) uui=x1i-x2i; void mtx_sub(double aa1d

11、imdim,double aa2dimdim,double mmdimdim) / 矩阵-矩阵 long i=0,j=0; for(i=0;idim;i+) for(j=0;jdim;j+) mmij=aa1ij-aa2ij; void mtx_mtx(double aa1dimdim,double aa2dimdim,double mmdimdim) / 矩阵* 矩阵 long i=0,j=0,k=0,t=0; for(i=0;idim;i+) for(j=0;jdim;j+) mmij=0; for(t=0;tdim;t+) mmij=mmij+aa1it*aa2tj; void vtr_

12、con(double xdim,double c,double uudim) / 向量乘、除 常数 long i=0; for(i=0;idim;i+) uui=xi/c; void mtx_con(double aadimdim,double c,double mmdimdim) / 矩阵 乘、除常数 long i=0,j=0; for(i=0;idim;i+) for(j=0;jdim;j+) mmij=mmij*c; void a_at(double aadimdim,double aatdimdim) / 矩阵转置 long i=0,j=0; for(i=0;idim;i+) for(

13、j=0;jdim;j+) aatij=aaji; void nssjh(double Ardimdim) / 拟上 三角化函数 double A1dimdim=0,0; double Hrdimdim=0,0; double srdim=0; double urdim=0; double prdim=0; double vrdim=0; double qrdim=0; double wrdim=0; double cr=0; double hr=0; double tr=0; long r=0,i=0,j=0; long i0=0; /第 r 列下半 部 0 的个数 for(r=0;rdim-

14、2;r+) for(i=r+2;idim;i+) / 判断 Ar 下半部分是否全为 0 if(Arir=0) i0+; else break; if(i0=dim-r-2) /若 Ar 下半部分全为 0，则 Hr 为 I for(i=0;idim;i+) for(j=0;jdim;j+) if(i=j) Hrij=1; else Hrij=0; else / 若 Ar 下 半部分不全为 0 ，对第 r 列拟上三角化 / 取 sr for(i=0;i=r;i+) sri=0; for(i=r+1;idim;i+) sri=Arir; / 取 cr if(Arr+1r=0) cr=length(s

15、r); else cr=0-length(sr); /计算 ur=sr-cr*er srr+1=srr+1-cr; /计算 hr hr=length(sr); hr=hr*hr/2; /计算 vr vtr_con(sr,hr,vr); /计算 qr mtx_vtr(Ar,vr,qr); /计算 pr a_at(Ar,A1); mtx_vtr(A1,sr,pr); /计算 tr tr=vtr_vtr(pr,vr); /计算 wr vtr_con(vr,1/tr,wr); vtr_sub(qr,wr,wr); /计算 Hr vtr_vtr2(sr,vr,Hr); mtx_sub(In,Hr,Hr)

16、; /计算 Q mtx_mtx(Q,Hr,Qf); for(i=0;idim;i+) for(j=0;jdim;j+) Qij=Qfij; /计算 A （r+1 ） vtr_vtr2(wr,sr,A1); mtx_sub(Ar,A1,Ar); vtr_vtr2(vr,pr,A1); mtx_sub(Ar,A1,Ar); i0=0; void qrfj(double Ardimdim) /QR 分 解函数 double A1dimdim=0,0; double Qrdimdim=0,0; double Rrdimdim=0,0; double Hrdimdim=0,0; double srdim

17、=0; double urdim=0; double prdim=0; double vrdim=0; double qrdim=0; double cr=0; double hr=0; long r=0,i=0,j=0; long i0=0; / 第 r 列 下半部 0 的个 数 for(i=0;idim;i+) for(j=0;jdim;j+) if(i=j) Qrij=1; else Qrij=0; for(r=0;rdim-1;r+) for(i=r+1;idim;i+) / 判断 Ar 下半部分是否全为 0 if(Arir=0) i0+; else break; if(i0=dim-

18、r-1) /若 Ar 下半部分全为 0，则 Hr 为 I for(i=0;idim;i+) for(j=0;jdim;j+) if(i=j) Hrij=1; else Hrij=0; else / 若 Ar 下 半部分不全为 0 ，对第 r 列拟上三角化 / 取 sr for(i=0;ir;i+) sri=0; for(i=r;idim;i+) sri=Arir; / 取 cr if(Arrr=0) cr=length(sr); else cr=0-length(sr); /计算 ur=sr-cr*er srr=srr-cr; /计算 hr hr=length(sr); hr=hr*hr/2;

19、 /计算 vr vtr_con(sr,hr,vr); /计算 qr mtx_vtr(Qr,sr,qr); /计算 pr a_at(Ar,A1); mtx_vtr(A1,sr,pr); /计算 Hr vtr_vtr2(sr,vr,Hr); mtx_sub(In,Hr,Hr); /计算 Qr vtr_vtr2(qr,vr,A1); mtx_sub(Qr,A1,Qr); /计算 Ar vtr_vtr2(vr,pr,A1); mtx_sub(Ar,A1,Ar); i0=0; for(i=0;idim;i+) /QR 分解 后把 Q 暂时 存入全局变量 FF for(j=0;jdim;j+) FFij=

20、Qrij; / 计算 Q mtx_mtx(Q,FF,Qf); for(i=0;idim;i+) for(j=0;jdim;j+) Qij=Qfij; /A(k+1)=R(k)*Q(k) mtx_mtx(Ar,Qr,A1); for(i=0;iA(k+1) mtx_mtx(Ak,FF,A1); a_at(FF,A2); mtx_mtx(A2,A1,Ak); void qr_sbwytz(double Ardimdim) /QR 分 解法求特征根 long k=0,i=0,j=0; long state=0; / 状态 long m=dim; double tzdim2=0,0; double b

21、=0,c=0; double aa; nssjh(Ar); while(m0) if(m2) if(fabs(Arm-1m-2)=er0) / 得到一个实 特征根 state=1; else if(fabs(Arm-2m-3)=er0) / 得到一对 共轭特征根 state=2; else if(k=L) / 判断是否 到最大次数 printf(“k=%dn“,k); m=0; else qr_sbwy(Ar); k=k+1; else if(m=1) state=1; else if(m=2) state=2; else if(m=0) printf(“failednk=%dnn“,k);

22、switch(state) case 1: / 得 到一个实特征根 tzi0=Arm-1m-1; printf(“tz%d=%12.12en“,dim-i,tzi0); i+; m-; state=0; break; case 2: /解一对共轭特征根 b=0-Arm-1m-1-Arm-2m-2; c=Arm-1m-1*Arm-2m-2-Arm-1m-2*Arm-2m-1; tzi0=-b/2; tzi1=-sqrt(4*c-b*b)/2; printf(“tz%d=%12.12e+%12.12en“,dim-i,tzi0,tzi1); i+; tzi0=-b/2; tzi1=sqrt(4*c

23、-b*b)/2; printf(“tz%d=%12.12e-%12.12en“,dim-i,tzi0,tzi1); i+; m=m-2; state=0; break; void init() long i=0,j=0; / 给出初 始矩阵 for(i=1;i=dim;i+) for(j=1;j=dim;j+) if(i=j) Ai-1j-1=1.52*cos(i+1.2*j); else Ai-1j-1=sin(0.5*i+0.2*j); / 单位阵 for(i=0;idim;i+) for(j=0;jdim;j+) if(i=j) Inij=1; else Inij=0; /Q 初始 值为

24、单位阵 for(i=0;idim;i+) for(j=0;jdim;j+) if(i=j) Qij=1; else Qij=0; （1 ） 拟上三角化 void main(void) init(); nssjh(A); printf(“A(n-1)=n“); display_m(A); （2 ） 求特征值和特征向量 void main(void) init(); qr_sbwytz(A); printf(“A(k+1)=n“); display_m(A); printf(“Q=n“); display_m(Q); 三、运行结果 （1 ） 拟上三角化 A(n-1)= -8.94521698228

25、1e-001 -9.933136491826e-002 -1.099831758877e+000 -7.66503 8709077e-001 1.707601141456e-001 -1.934882558889e+000 -8.390208705245e -002 9.132565113143e-001 -6.407977009188e-001 1.946733678685e-001 -2.347878362416e+000 2.372057921598e+000 1.827998552316e+000 3.266556 884714e-001 2.082360583635e-001 2.0

26、88987009941e+000 1.847861910289e- 001 -1.263015266080e+000 6.790694668499e-001 -4.672150886500e-001 -3.400178888237e-016 1.735954469946e+000 -1.165023367477e+000 -1.24674 4443518e+000 -6.298225489084e-001 -1.984820180992e+000 2.975750060800e- 001 6.339300596595e-001 -1.308518928772e-001 3.0403010360

27、95e-001 -1.169346032160e-016 1.952009050168e-017 -1.292937563924e+000 -1.12623 9225902e+000 1.190782911924e+000 -1.308772983895e+000 1.860151662666e- 001 4.236733936881e-001 -1.019600826545e-001 1.943660914505e-001 -1.603077371689e-017 5.393945153379e-017 -1.533762873820e-016 1.577711 153032e+000 8.

28、169358328160e-001 4.461531723827e-001 -4.365092541609e -002 -4.665979167188e-001 2.941231566184e-001 -1.034421113665e-001 3.131441504873e-016 3.600238995442e-017 9.033164589718e-017 9.171310 651350e-018 -7.728975134989e-001 -1.601028244046e+000 -2.912685474826e -001 -2.434337858321e-001 6.7362860845

29、10e-001 2.624772904937e-001 3.016667459975e-016 -2.658831366350e-016 -6.526288055706e-017 -8.80991 5803872e-017 3.378335595382e-017 -7.296773946362e-001 -7.965456279828e -003 9.710739102007e-001 -1.298967368574e-001 2.780242081241e-002 -2.559330035627e-016 -1.244642736267e-016 1.682913291677e-016 -5

30、.02328 8971095e-018 4.125975483652e-017 -5.536518281074e-017 7.945539612976e- 001 -4.525143454606e-001 5.048901527575e-001 -1.211210193512e-001 6.554311307585e-017 1.143144843268e-016 -1.341950015062e-016 -6.65980 9074964e-017 -2.763275502874e-017 -1.001176808787e-018 0.000000000000e+ 000 7.03991137

31、3514e-001 1.267535523498e-001 -3.714696735513e-001 7.887019170884e-017 -1.824459988867e-016 -4.119393538865e-017 -2.23825 3293760e-016 1.602015267643e-017 -3.896419940204e-018 0.000000000000e+ 000 0.000000000000e+000 -4.919586872214e-001 4.081509766399e-001 （2 ） 求特征值和特征向量 tz10=4.954990923632e-002 tz

32、9=6.489488202108e-001 tz8=9.432879572769e-001 tz7=-9.891143464723e-001+-1.084758631501e-001 tz6=-9.891143464723e-001-1.084758631501e-001 tz5=-1.493147080915e+000 tz4=1.590313458807e+000 tz3=-2.336865932239e+000+-8.934379210213e-001 tz2=-2.336865932239e+000-8.934379210213e-001 tz1=3.389613438816e+000

33、 A(k+1)= 3.389613438816e+000 -6.760784761554e-001 1.061018160428e+000 8.837866 591599e-002 2.662196283780e-001 1.511975111218e+000 1.176722809109e+ 000 -1.076059260864e-001 1.004709088549e+000 4.947584627475e-001 7.172433346708e-072 -2.576495582619e+000 -2.344292534990e+000 -7.99570 4025565e-002 4.6

34、39106523634e-002 -2.221472389689e+000 -3.992548977849e -001 -1.790363956051e-001 -1.862039486392e+000 -9.550111925940e-001 3.108444118593e-086 3.649944174157e-001 -2.097236281858e+000 6.347483 410359e-001 -1.292694251658e-002 -1.314495195394e+000 1.111156715258e- 001 -2.862368601443e-001 -9.09778366

35、6310e-001 -8.144941738974e-002 1.609609990678e-193 9.445945194823e-139 1.541740815009e-107 1.590313 458807e+000 -1.546422732140e-002 -6.884638866309e-001 -2.257044474813e -001 -2.691095898199e-002 -5.129587352611e-001 -2.744100394559e-001 -8.893280172785e-208 -2.418435398724e-137 1.558547256956e-136

36、 2.930360 739935e-015 -1.493147080915e+000 -9.768514017102e-002 -5.050320228930e -002 6.389148148084e-003 -3.476259059528e-002 -4.492865707191e-002 2.232104685347e-303 1.389609222842e-232 9.133555102531e-233 1.251162 431503e-125 -7.360755271166e-097 -7.356906393162e-001 5.863881195154e- 001 -2.63026

37、9467690e-001 1.138644581274e-001 3.864695938656e-001 -3.764772518155e-304 -4.745036630717e-233 -3.061768285351e-233 -3.58620 9717329e-126 1.675549200756e-112 -1.295909410609e-001 -1.242538053628e +000 -9.191699255585e-002 -4.225065931832e-001 -2.227969665138e-001 2.209399562964e-316 -3.312336222183e

38、-246 -3.212420877712e-246 -1.17654 9574627e-138 1.619095130415e-123 8.330371159134e-028 8.162965389350e- 013 9.432879572769e-001 1.889770633289e-001 1.403752464446e-001 -1.581010066692e-322 5.928787750095e-323 -1.482196937524e-322 -2.03587 5909915e-227 1.160837712045e-212 2.399148143087e-116 3.410359343893e- 116 -5.683390400115e-090 6.48