北航数值分析报告大作业二90761
- 格式：doc
- 大小：1.14 MB
- 文档页数：17
数值分析第二题1 算法设计方案要想得出该题的答案首先要将矩阵A 进行拟上三角化，把矩阵A 进行QR 分解。
要得出矩阵A 的全部特征值先对A 进行QR 的双步位移得出特征值。
采用列主元的高斯消元法求解特征向量。
1.1 A 的拟上三角化因为对矩阵进行QR 分解并不改变矩阵的结构，因此在进行QR 分解前对矩阵A 进行拟上三角化可以大大减少计算机的计算量，提高程序的运行效率。
具体算法如下所示，记A A =)1(，并记)(r A 的第r 列至第n 列的元素为()n r r j n i a r ij,,1,;,,2,1)( +==。
对于2,,2,1-=n r 执行1. 若()n r r i a r ir,,3,2)( ++=全为零，则令)()1(r r A A =+,转5；否则转2。
2. 计算()∑+==nr i r irr a d 12)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若 )(,12r rr r r r a c c h +-=3. 令()n Tr nrr r r r r r r r R a a c a u ∈-=++)()(,2)(,1,,,,0,,0 。
4. 计算r r T r r h u A p /)(= r r r r h u A q /)(=r r Tr r h u p t /=r r r r u t q -=ωTrr T r r r r p u u A A --=+ω)()1(5. 继续。
1.2 A 的QR 分解具体算法如下所示，记)1(1-=n A A ，并记[]nn r ij r a A ⨯=)(，令I Q =1 对于1,,2,1-=n r 执行1.若()n r r i a r ir,,3,1)( ++=全为零，则令r r Q Q =+1r r A A =+1,转5；否则转2。
2.计算()∑==nri r irr a d 2)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若)(,2r r r r r r a c c h -=3.令()n Tr nrr r r r r r r r R a a c a u ∈-=+)()(,2)(,,,,,0,,0 。
4.计算r r T r r h u A p /)(= r r r r h u A q /)(=r r Tr r h u p t /=r r r r u t q -=ωT rr T r r r r p u u A A --=+ω)()1(5.继续。
1. 3 A 的特征值为了加快收敛速度，采用带双步移位的QR 分解求解A 的全部特征值，具体算法如下，(1) 使用矩阵的拟上三角化算法把矩阵A 化为你上三角阵；给定经度水平和最大迭代次数L 。
(2) 记，令k=1,m=n 。
(3) 如果，则得到A 的一个特征值，置m:=m-1(降阶)，转(5)，否则转(4)。
(4) 如果，则得到矩阵A 的两个特征值,为二阶子阵⎥⎥⎦⎤⎢⎢⎣⎡=----)()(1,(k),1)(1,1 k mm k m m m m k m m k a a a a D 的两个特征根，置m:=m-1(降阶)，转(5)；否则转(6)。
(5) 如果m=1，则得到A 的一个特征根，转(9)，如果m=0，转(9)。
否则转(3)。
(6) 如果k=L ，则计算失败。
否则转(7)。
(7) 记，计算kk T k k k k k k k k k k mm k m m k mm k m m k mmk m m Q A Q A QR M R Q M I tI sA A M a a a a t a a s ==+-=-=+=+------12)(,1)(1,)()(1,1)()(1,1)分解作对阶单位矩阵是( )m ( (8) ，转(3)。
(9) A 的全部特征值以计算完毕，停止计算。
其中，k M 的QR 分解与1+k A 的计算用下列算法实现：记k m m r ij r m m ij k A C b B b M B ====⨯⨯1)()1(1,][,][。
对于1,,2,1-=m r 执行1. 若()m r r i b r ir,,2,1)( ++=全为零，则令r r r r C C B B ==++11,,转5；否则转2。
2. 计算()∑==mri r irr b d 2)(()()r r r r r r r rr r d c a d b c ==-=+则取若,0sgn )(,1)( )(2r rr r r r b c c h -=3. 令()m Tr mrr r r r r rr r R b b c b u ∈-=+)()(,1)(,,,,0,,0 。
4. 计算r r Tr r h u B v /=T r r r r v u B B -=+1r r T r r h u C p /=r r r r h u C q /= r r T r r h u p t /=r r r r u t q -=ωT r r T r r r r p u u C C --=+ω15. 继续。
此算法执行完后，就得到m k C A =+1。
1.4 A 的特征向量记λ为矩阵A 的实特征值，x 为对应的特征向量。
则Ax=λx ，即(A-λI)x=0的解即为矩阵A 的特征向量。
因此对于特征向量的求解可采用列主元素的Gauss 消元法。
其具体算法如下，1. 消元过程对于1,,2,1-=n k 执行(1) 选行号k i ，使)()(max k ikni k k ki a a k ≤≤=。
(2) 交换)(k kj a 与),,1,()(n k k j a k ji k +=所含的数值。
(3) 对于n k k i ,,2,1 ++=计算),,2,1( /)()()1()()(n k k j a m a a a a m k ki ik k ij k ijk kkk ik ik ++=-==+2. 回代过程()1,,2,1 1)(1)( --=⎪⎪⎭⎫ ⎝⎛-==∑+=n n k a x a x x k kk nk j j k kj kn最终得到的向量T n x x x )1,,,(21= 的即为A 对应于实特征值λ的特征向量。
2.程序运行结果2.1拟上三角化后2.2 QR分解Q阵R阵RQ阵2.3QR的双步位移求特征值和特征向量特征值特征向量3源程序#include<stdio.h>#include<math.h>#include<stdlib.h>#define L 2500#define n 11#define E 1e-12int i,j,s,p,k,ik,nR,nC;doubleA[n][n],q[n][n],r[n][n],rq[n][n],I[n][n],tzR[n],tzC[2][n],M[n][n],v[n ];double P[n],W[n],u[n],Q[n];double dr,cr,hr,ar,tr,s1,t,x,sum;void createA()////生成矩阵A{for(i=1;i<11;i++){for(j=1;j<11;j++){if(j!=i)A[i][j]=sin(0.5*i+0.2*j);elseA[i][j]=1.5*cos(i+1.2*j);}}}void nszj()//拟上三角化A{for(s=1;s<n-2;s++){for(ar=0.0,i=s+2;i<n;i++)ar+=A[i][s]*A[i][s];if(ar==0)continue;else {ar+=A[s+1][s]*A[s+1][s];dr=sqrt(ar);if(A[s+1][s]>0) cr=-dr;else cr=dr;hr=cr*cr-cr*A[s+1][s];for(i=1;i<=s;i++)u[i]=0.0;u[s+1]=A[s+1][s]-cr;for(i=s+2;i<n;i++)u[i]=A[i][s];for(j=1;j<n;j++){for(P[j]=0.0,i=1;i<n;i++) P[j]+=A[i][j]*u[i]/hr;}for(tr=0.0,i=1;i<n;i++){tr+=P[i]*u[i]/hr;}for(i=1;i<n;i++){for(Q[i]=0.0,j=1;j<n;j++)Q[i]+=A[i][j]*u[j]/hr;}for(i=1;i<n;i++){W[i]=Q[i]-tr*u[i];}for(i=1;i<n;i++){for(j=1;j<n;j++)A[i][j]-=(W[i]*u[j]+u[i]*P[j]);}}}for(i=1;i<n;i++){for(j=1;j<n;j++)if(fabs(A[i][j])<E)A[i][j]=0.0;}printf("拟上三角化后A(n-1):\n");for(i=1;i<n;i++){for(j=1;j<n;j++)printf("%1.12e,",A[i][j]);printf("\n");}}void QR()//QR分解{double u[n],w[n],F[n];for(s=1;s<n;s++)for(p=1;p<n;p++)r[s][p]=A[s][p];for(s=1;s<n-1;s++){for(dr=0.0,i=s;i<n;i++)dr+=r[i][s]*r[i][s];dr=sqrt(dr);if (dr==0)continue;else{if(A[s][s]>0) cr=-dr;else cr=dr;hr=cr*cr-cr*r[s][s];for(i=1;i<s;i++)u[i]=0;u[s]=r[s][s]-cr;for(i=s+1;i<n;i++)u[i]=r[i][s];for(i=1;i<n;i++){for (F[i]=0.0,j=1;j<n;j++)F[i]+=r[j][i]*u[j]/hr;}for(i=1;i<n;i++){for(j=1;j<n;j++)r[i][j]=r[i][j]-u[i]*F[j];}for(i=1;i<n;i++){for(w[i]=0.0,j=1;j<n;j++)w[i]+=q[i][j]*u[j];}for(i=1;i<n;i++){for(j=1;j<n;j++)q[i][j]-=w[i]*u[j]/hr;}}}for(i=1;i<n;i++){for(j=1;j<n;j++){for(rq[i][j]=0.0,s=1;s<n;s++)rq[i][j]+=r[i][s]*q[s][j];}}printf("生成的Q阵：\n");for(i=1;i<n;i++){for(j=1;j<n;j++)if(fabs(q[i][j])<E)q[i][j]=0.0;}for(i=1;i<n;i++){for(j=1;j<n;j++)printf("%1.12e,",q[i][j]);printf("\n");}printf("生成的R阵：\n");for(i=1;i<n;i++){for(j=1;j<n;j++)if(fabs(r[i][j])<E)r[i][j]=0.0;}for(i=1;i<n;i++){for(j=1;j<n;j++)printf("%1.12e,",r[i][j]);printf("\n");}printf("生成的RQ阵：\n");for(i=1;i<n;i++){for(j=1;j<n;j++)if(fabs(rq[i][j])<E)rq[i][j]=0.0;}for(i=1;i<n;i++){for(j=1;j<n;j++)printf("%1.12e,",rq[i][j]);printf("\n");}}void QRtwometh()//QR双步位移求特征值{int K=1,m=10;nR=0,nC=0;loop3:if(fabs(A[m][m-1])<=E){nR++;tzR[nR]=A[m][m];m--;goto loop4;}else goto loop5;loop4:if(m==1){nR++;tzR[nR]=A[1][1];goto loop9;}else if(m==2){s1=A[m-1][m-1]+A[m][m];t=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m];x=s1*s1-4*t;if(x>=0){nR++;tzR[nR]=(s1+sqrt(x))/2;nR++;tzR[nR]=(s1-sqrt(x))/2;}else{nC++;tzC[0][nC]=s1/2;tzC[1][nC]=sqrt(-x)/2;nC++;tzC[0][nC]=s1/2;tzC[1][nC]=-sqrt(-x)/2;}goto loop9;}else goto loop3;loop5:{if(fabs(A[m-1][m-2])<=E){s1=A[m-1][m-1]+A[m][m];t=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m];x=s1*s1-4*t;if(x>=0){nR++;tzR[nR]=(s1+sqrt(x))/2;nR++;tzR[nR]=(s1-sqrt(x))/2;}else{nC++;tzC[0][nC]=s1/2;tzC[1][nC]=sqrt(-x)/2;nC++;tzC[0][nC]=s1/2;tzC[1][nC]=-sqrt(-x)/2;}m--;m--;goto loop4;}else goto loop6;}loop6:{if(K==L)printf("计算终止，未能得到全部特征值\n");else goto loop7;}loop7:{s1=A[m-1][m-1]+A[m][m];t=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m];for(i=1;i<=m;i++)for(j=1;j<=m;j++){for(sum=0.0,s=1;s<=m;s++)sum+=A[i][s]*A[s][j];M[i][j]=sum-s1*A[i][j]+t*I[i][j];}for(s=1;s<=m;s++){for(ar=0.0,i=s+1;i<=m;i++)ar+=M[i][s]*M[i][s];if(ar==0)continue;else{ar+=M[s][s]*M[s][s];dr=sqrt(ar);if(M[s][s]>0)cr=-dr;else cr=dr;hr=cr*cr-cr*M[s][s];for(i=1;i<s;i++)u[i]=0.0;u[s]=M[s][s]-cr;for(i=s+1;i<=m;i++)u[i]=M[i][s];for(j=1;j<=m;j++)for(v[j]=0.0,i=1;i<=m;i++)v[j]+=M[i][j]*u[i]/hr;for(i=1;i<=m;i++)for(j=1;j<=m;j++)M[i][j]-=u[i]*v[j];for(j=1;j<=m;j++)for(P[j]=0.0,i=1;i<=m;i++)P[j]+=A[i][j]*u[i]/hr;for(i=1;i<=m;i++)for(Q[i]=0.0,j=1;j<=m;j++)Q[i]+=A[i][j]*u[j]/hr;for(tr=0.0,i=1;i<=m;i++)tr+=P[i]*u[i]/hr;for(i=1;i<=m;i++)W[i]=Q[i]-tr*u[i];for(i=1;i<n;i++)for(j=1;j<n;j++)A[i][j]-=(W[i]*u[j]+u[i]*P[j]);}}goto loop8;}loop8:{k++;goto loop3;}loop9: ;printf("矩阵的全部特征值为：\n");for(i=1;i<=nR;i++)printf("%1.12e\n",tzR[i]);for(i=1;i<=nC;i++){printf("%1.12e ",tzC[0][i]);if(tzC[1][i]>=0)printf("%1.12e\n",tzC[1][i]);else printf("%1.12e\n",tzC[1][i]);}}void gas()//列主元的高斯消元法求解特征向量{double ch,m[n],x[n][n];for(p=1;p<=nR;p++){createA();for(i=1;i<n;i++)for(j=1;j<n;j++)A[i][j]-=tzR[p]*I[i][j];for(k=1;k<n;k++){for(ik=k,i=k;i<n;i++)if(A[ik][k]<A[i][k])ik=i;for(j=k;j<n;j++){ch=A[ik][j];A[ik][j]=A[k][j];A[k][j]=ch;}for(i=k+1;i<n;i++){m[i]=A[i][k]/A[k][k];for(j=k+1;j<n;j++)A[i][j]-=m[i]*A[k][j];}}x[p][n-1]=1.0;for(k=n-2;k>0;k--){for(ar=0.0,j=k+1;j<=n;j++)ar+=A[k][j]*x[p][j];}}for(p=1;p<=nR;p++){printf("对应特征值%1.12e的特征向量为:\n",tzR[p]);for(i=1;i<n;i++)printf("%1.12e\n",x[p][i]);}}void main(){for(i=1;i<n;i++){for(j=1;j<n;j++){if (i==j){I[i][j]=1;q[i][j]=1;}else{I[i][j]=0;q[i][j]=0;}}}createA();nszj();QR();QRtwometh();gas();}。
