北航数值分析大作业二.doc

数值分析大作业（二）学院名称宇航学院专业名称航空宇航推进理论与工程学生姓名段毓学号SY16153062016年11月5日1 算法设计方案首先将矩阵A进行拟上三角化，把矩阵A进行QR分解，计算出RQ。要得出矩阵A的全部特征值，首先对A进行QR的双步位移得出特征值。最后，采用列主元的高斯消元法求解特征向量。1.1 A的拟上三角化因为对矩阵进行QR分解并不改变矩阵的结构，因此在进行QR分解前对矩阵A进行拟上三角化可以大大减少计算机的计算量，提高程序的运行效率。具体算法如下所示，记，并记的第r列至第n列的元素为。对于执行若全为零，则令,转5；否则转2。计算令。计算 继续。1.2 A的QR分解 具体算法如下所示，记，并记，令对于执行1.若全为零，则令,转5；否则转2。2.计算3.令。4.计算 5.继续。1. 3 A的特征值为了加快收敛速度，采用带双步移位的QR分解求解A的全部特征值，具体算法如下，使用矩阵的拟上三角化算法把矩阵A化为你上三角阵；给定经度水平和最大迭代次数L。记，令k=1,m=n。如果，则得到A的一个特征值，置m:=m-1(降阶)，转(5)，否则转(4)。如果，则得到矩阵A的两个特征值,为二阶子阵的两个特征根，置m:=m-1(降阶)，转(5)；否则转(6)。如果m=1，则得到A的一个特征根，转(9)，如果m=0，转(9)。否则转(3)。如果k=L，则计算失败。否则转(7)。记，计算，转(3)。A的全部特征值以计算完毕，停止计算。其中，的分解与的计算用下列算法实现：记。对于执行若全为零，则令,转5；否则转2。计算令。计算继续。此算法执行完后，就得到。1.4 A的特征向量记为矩阵A的实特征值，x为对应的特征向量。则Ax=x，即(A-I)x=0的解即为矩阵A的特征向量。因此对于特征向量的求解可采用列主元素的Gauss消元法。其具体算法如下，消元过程对于执行选行号，使。交换与所含的数值。对于计算回代过程最终得到的向量的即为对应于实特征值的特征向量。2. 程序运行结果2.1拟上三角化后2.2 QR分解Q阵R阵RQ阵2.3QR的双步位移求特征值和特征向量特征值特征向量3源程序#include stdio.h#include math.h#define n 11#define E 1e-12#define L 2500double fuzhi(double a1010);/给矩阵A赋值double nishangsanjiaohua(double a1010);double QR(double a1010);double sgn(double x);void gauss();double A1010;double tzRn,Inn;int nR;void QRshuangbu();void main()fuzhi(A);nishangsanjiaohua(A);QR(A);QRshuangbu();gauss(); double fuzhi(double a1010)/给矩阵A赋值 int i, j; for(i = 1; i 11; i+) for(j = 1; j 11; 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); double sgn(double x)/定义sgn函数 if(x 0)return 1;if(x = 0)return -1;double nishangsanjiaohua(double a1010)int r, i, j;double pingfanghe;double drpingfang;double dr, cr, hr;double ur10, pr10, qr10, tr, wr10;for(r = 1; r 9; r+)pingfanghe = 0;drpingfang = 0;for(i = 1; i 11; i+)uri-1 = 0;pri-1 = 0;qri-1 = 0;wri-1 = 0;tr = 0;/置零 for(i = r + 2; i 11; i+)pingfanghe = pingfanghe + ai-1r-1 * ai-1r-1;if(pingfanghe = 0)continue;/(1)elsefor(i = r + 1; i 11; i+)drpingfang = drpingfang + ai-1r-1 * ai-1r-1;dr = sqrt(drpingfang);/(2.1)drcr = -sgn(arr-1) * dr;/(2.2)crhr = cr * cr - cr * arr-1;/(2.3)hrfor(i = r + 1; i 11; i+)uri-1 = ai-1r-1;urr = urr - cr;/(3)urfor(j = 1; j 11; j+)for(i = 1; i 11; i+)prj-1 = prj-1 + ai-1j-1 * uri-1;prj-1 = prj-1 / hr;/(4.1)prfor(i = 1; i 11; i+)for(j = 1; j 11; j+)qri-1 = qri-1 + ai-1j-1 * urj-1;qri-1 = qri-1 / hr;/(4.2)qrfor(i = 1; i 11; i+)tr = tr + pri-1 * uri-1;tr = tr / hr;for(i = 1; i 11; i+)wri-1 = qri-1 - tr * uri-1;for(i = 1; i 11; i+)for(j = 1; j 11; j+)ai-1j-1 = ai-1j-1 - wri-1 * urj-1 - uri-1 * prj-1;for(i = 1; i 11; i+)for(j = 1; j 11; j+)printf(A%d,%d=%.12et, i, j, ai-1j-1);printf(n);double QR(double a1010)double Q1010;double RQ1010;int i, j;int t;int r;double pingfanghe;double drpingfang;double dr, cr, hr;double ur10, wr10, pr10;for(i = 1; i 11; i+)for(j = 1; j 11; j+)if(i = j)Qi-1j-1 = 1;elseQi-1j-1 = 0;for(r = 1; r 10; r+)pingfanghe = 0;drpingfang = 0;for(i = 1; i 11; i+)uri-1 = 0;wri-1 = 0;pri-1 = 0;for(i = r + 1; i 11; i+)pingfanghe = pingfanghe + ai-1r-1 * ai-1r-1;if(pingfanghe = 0)continue;elsefor(i = r; i 11; i+)drpingfang = drpingfang + ai-1r-1 * ai-1r-1;dr = sqrt(drpingfang);cr = - sgn(ar-1r-1) * dr;hr = cr * cr - cr * ar-1r-1;for(i = r; i 11; i+)uri-1 = ai-1r-1;urr-1 = urr-1 - cr;for(i = 1; i 11; i+)for(j = 1; j 11; j+)wri-1 = wri-1 + Qi-1j-1 * urj-1;for(i = 1; i 11; i+)for(j = 1; j 11; j+)Qi-1j-1 = Qi-1j-1 - wri-1 * urj-1 / hr;for(i = 1; i 11; i+)for(j = 1; j 11; j+)pri-1 = pri-1 + aj-1i-1 * urj-1 / hr;for(i = 1; i 11; i+)for(j = 1; j 11; j+)ai-1j-1 = ai-1j-1 - uri-1 * prj-1;printf(正交矩阵Q：n);for(i = 1; i 11; i+)for(j = 1; j 11; j+)printf(%.12et,Qi-1j-1);printf(n); printf(上三角矩阵R:n);for(i = 1; i 11; i+)for(j = 1; j 11; j+)printf(%.12et,ai-1j-1);printf(n);for(i = 1; i 11; i+)for(j = 1; j 11; j+)for(t = 1; t 11; t+)RQi-1j-1 = RQi-1j-1 + ai-1t-1 * Qt-1j-1;/计算RQ printf(RQ阵：n);for(i = 1; i 11; i+)for(j = 1; j 11; j+)printf(%.12et, RQi-1j-1);printf(n);void QRshuangbu()/QR双步位移求特征值int K=1,m=10;int i,j,s,p,k,ik,nC;double qnn,rnn,rqnn,tzC2n,Mnn,vn;double Pn,Wn,un,Qn;double dr,cr,hr,ar,tr,s1,t,x,sum;nR=0,nC=0;loop3:if(fabs(Amm-1)=E) nR+; tzRnR=Amm; m-; goto loop4; else goto loop5;loop4:if(m=1) nR+; tzRnR=A11; goto loop9; else if(m=2) s1=Am-1m-1+Amm; t=Am-1m-1*Amm-Amm-1*Am-1m; x=s1*s1-4*t; if(x=0) nR+; tzRnR=(s1+sqrt(x)/2; nR+; tzRnR=(s1-sqrt(x)/2; else nC+; tzC0nC=s1/2; tzC1nC=sqrt(-x)/2; nC+; tzC0nC=s1/2; tzC1nC=-sqrt(-x)/2; goto loop9; else goto loop3;loop5: if(fabs(Am-1m-2)=E) s1=Am-1m-1+Amm; t=Am-1m-1*Amm-Amm-1*Am-1m; x=s1*s1-4*t; if(x=0) nR+; tzRnR=(s1+sqrt(x)/2; nR+; tzRnR=(s1-sqrt(x)/2; else nC+; tzC0nC=s1/2; tzC1nC=sqrt(-x)/2; nC+; tzC0nC=s1/2; tzC1nC=-sqrt(-x)/2; m-; m-; goto loop4; else goto loop6; loop6:if(K=L) printf(计算终止，未能得到全部特征值n); else goto loop7; loop7: s1=Am-1m-1+Amm; t=Am-1m-1*Amm-Amm-1*Am-1m; for(i=1;i=m;i+) for(j=1;j=m;j+) for(sum=0.0,s=1;s=m;s+) sum+=Ais*Asj; Mij=sum-s1*Aij+t*Iij; for(s=1;s=m;s+) for(ar=0.0,i=s+1;i=m;i+) ar+=Mis*Mis; if(ar=0) continue; else ar+=Mss*Mss; dr=sqrt(ar); if(Mss0) cr=-dr; else cr=dr; hr=cr*cr-cr*Mss; for(i=1;is;i+) ui=0.0; us=Mss-cr; for(i=s+1;i=m;i+) ui=Mis; for(j=1;j=m;j+) for(vj=0.0,i=1;i=m;i+) vj+=Mij*ui/hr; for(i=1;i=m;i+) for(j=1;j=m;j+) Mij-=ui*vj; for(j=1;j=m;j+) for(Pj=0.0,i=1;i=m;i+) Pj+=Aij*ui/hr; for(i=1;i=m;i+) for(Qi=0.0,j=1;j=m;j+) Qi+=Aij*uj/hr; for(tr=0.0,i=1;i=m;i+) tr+=Pi*ui/hr; for(i=1;i=m;i+) Wi=Qi-tr*ui; for(i=1;in;i+) for(j=1;jn;j+) Aij-=(Wi*uj+ui*Pj); goto loop8; loop8: k+; goto loop3; loop9: ;printf(矩阵的全部特征值为：n); for(i=1;i=nR;i+) printf(%1.12en,tzRi); for(i=1;i=nC;i+) printf(%1.12e ,tzC0i); if(tzC1i=0) printf(%1.12en,tzC1i); else printf(%1.12en,tzC1i); void gauss()int N = 10;double BN;int k, i, j;int s, t, dijige;double m;double temp;int ik;double store;for(dijige = 1; dijige 11; dijige+)for(i = i; i 11; i+)Bi-1 = 0;Ai-1i-1 = Ai-1i-1 - tzRdijige-1;for(k = 1; k N; k+)/选行号 temp = Ak-1k-1;for(i = k; i N + 1; i+)ik = k;if (temp abs(Ai-1k-1)temp = abs(Ai-1k-1);ik = i;/交换akj与aikj for(j = k; j N + 1; j+)store = Aik-1j-1;Aik-1j-1 = Ak-1j-1;Ak-1j-1 = store;/交换bk与bikstore = Bi