数值分析第二次大作业史立峰SY1505327一、 方案（1）利用循环结构将sin(0.50.2)()1.5cos( 1.2)(){i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值，得到需要变换的矩阵A ；（2）然后，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵A (n -1)。

对A 拟上三角化，得到拟上三角矩阵A (n -1)，具体算法如下： 记A(1)=A ，并记A(r)的第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)( ++=全为零，则令A(r+1) =A(r),转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. 令()nTr 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 -=ωT rr T r r r r p u u A A --=+ω)()1(5. 继续。

（3）使用带双步位移的QR 方法计算矩阵A (n -1)的全部特征值，也是A 的全部特征值，具体算法如下：1. 给定精度水平0>ε和迭代最大次数L 。

2. 记n n ij n a A A ⨯-==][)1()1()1(，令n m k ==,1。

3. 如果ε≤-)(1,k m m a ，则得到A 的一个特征值)(k mm a ，置1:-=m m （降阶），转4；否则转5。

4. 如果1=m ，则得到的一个特征值)(11k a ，转11；如果1>m ，则转3。

5. 求2阶子阵⎥⎥⎦⎤⎢⎢⎣⎡=----)()(1,(k),1)(1,1 k mm k m m m m k m m k a a a a D的两个特征值1s 和2s ，即计算二次方程det )()()(1.12=++---k k mm k m m D a a λλ的两个根1s 和2s 。

6. 如果2=m ，则得到A 的两个特征值1s 和2s ，转11；否则转7。

7. 如果ε≤--)(2,1k m m a ，则得到A 的两个特征值1s 和2s ，置2:-=m m （降阶），转4；否则转88. 如果L k =，则计算终止，未得到A 的全部特征值；否则转9。

9. 记),1(][)()(m j i a A m m k ij k ≤≤=⨯，计算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 ( 10. 置1:+=k k ，转3。

11. 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 rrr r r b c c h -=3. 令()m Tr mr r r r r r rr r R b b c b u ∈-=+)()(,1)(,,,,0,,0 。

4. 计算r r T r r h u B v /=Trr 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 Tr r h u p t /=r r r r u t q -=ωT rr T r r r r p u u C C --=+ω15. 继续。

此算法执行完后，就得到m k C A =+1。

（4）计算Q,R ，一边求R*Q 矩阵。

记阶单位阵）n （令,][，并记A 1)(r 1I Q a A A n n r ij ===⨯ 对于1n ,,2,1-= r 执行 1.若()n r r i a r ir,,3,2)( ++=全为零，则令rr r r A A Q Q ==++11转5；否则转 2.2.计算()∑==mri r irr b d 2)(()()r r r r r r r rrr d c a d b c ==-=+则取,0sgn )(,1)(若)(2r rrr r r b c c h -=3.令()m Tr mr r r r r r rr r R b b c b u ∈-=+)()(,1)(,,,,0,,0 。

4.计算Trr r r rrTr rT r r r r r r p u A A h u A p h u Q Q u Q -==-==++1r 1r ωω5.继续当此算法执行完毕后，就得到正交矩阵n Q Q =和上三角矩阵R 6.然后计算出Q R 矩阵（5）用列主元素Gauss 消去法计算矩阵A 对应于实特征值的特征向量，具体算法如下：记)( )1(的实特征值为A I A Aλλ-=1. 消元过程对于1,,2,1-=n k 执行 (1) 选行号k i ，使)()(max k ikni k k k i aa k ≤≤=。

(2) 交换)(k kj a 与),,1,()(n k k j a k j i 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 ij k kkk ik ik ++=-==+2. 回代过程()1,,2,1 1)(1)( --=⎪⎪⎭⎫ ⎝⎛-==∑+=n n k a x ax x k kk nk j j k kjk n最终得到的向量Tn x x x )1,,,(21= 的即为A 对应于实特征值λ的特征向量。

二、源程序#include<iostream.h>#include<math.h>#include<iomanip.h>#define n 10#define E 1.0e-12void Householder(double a[n][n]);//拟上三角化函数double sgn(double a);//符号函数void QRfenjie(double a[n][n],double Q[n][n],double R[n][n]);void QR(double a[n][n],double L[n][2]);//带双步位移的QR分解void MxM(double M[n][n],double A[n][n],double B[n][n],int m);//矩阵相乘void solve(double a[n][n],double s1[2],double s2[2],int m);//解方程函数void Gauss(double a[n][n], double x[n]);//定义列主元高斯消去法函数void main(){int i,j,k;double a[n][n],qr[n][n],q[n][n],r[n][n],L[n][2],x[n];for(i=0;i<n;i++)//录入要进行变换的矩a，并对q初始赋值。

{for(j=0;j<n;j++){if(i==j)a[i][j]=1.5*cos(i+1+1.2*(j+1));elsea[i][j]=sin(0.5*(i+1)+0.2*(j+1));}}// cout<<endl;Householder(a);//调用拟上三角化函数cout<<endl<<endl<<"对矩阵A拟上三角化前七列的结果:"<<endl;for(i=0;i<n;i++){//k=0;//为了显示方便，每行显示三个元素，使用k来实现//cout<<"矩阵A的第"<<i+1<<"行元素为："<<endl;for(j=0;j<7;j++){if (fabs(a[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<a[i][j]<<"";//k++;//cout<<endl;//if(k%3==0)//判断每行是否到达三个元素}cout<<endl;}cout<<"对矩阵A拟上三角化后三列的结果:"<<endl;for(i=0;i<n;i++){//k=0;//为了显示方便，每行显示三个元素，使用k来实现//cout<<"矩阵A的第"<<i+1<<"行元素为："<<endl;for(j=7;j<n;j++){if (fabs(a[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<a[i][j]<<"";//k++;//if(k%3==0)//判断每行是否到达三个元素//cout<<endl;}cout<<endl;}cout<<endl;QRfenjie(a,q,r);QR(a,L);cout<<endl<<"对矩阵A进行QR分解后所得矩阵前七列的结果:"<<endl; for(i=0;i<n;i++){for(j=0;j<7;j++){if (fabs(a[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<a[i][j]<<"";}cout<<endl;}cout<<"对矩阵A进行QR分解后所得矩阵三列的结果:"<<endl;for(i=0;i<n;i++){for(j=7;j<n;j++){if (fabs(a[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<a[i][j]<<"";}cout<<endl;}/*cout<<"对矩阵A进行QR分解后所得矩阵:"<<endl;for(i=0;i<n;i++){k=0;cout<<"矩阵A的第"<<i+1<<"行元素为："<<endl;for(j=0;j<n;j++){if (fabs(a[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<a[i][j]<<"";k++;if(k%3==0)cout<<endl;}cout<<endl;}*/for(i=0;i<n;i++){for(j=0;j<n;j++)for(k=0,qr[i][j]=0;k<n;k++)qr[i][j]=qr[i][j]+r[i][k]*q[k][j];}cout<<endl<<"R*Q相乘的前七列:"<<endl;for(i=0;i<n;i++){//k=0;//cout<<"R*Q的第"<<i+1<<"行元素为："<<endl;for(j=0;j<7;j++){if (fabs(qr[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<qr[i][j]<<"";//k++;//if(k%3==0)//cout<<endl;}cout<<endl;}cout<<endl<<"R*Q相乘的后三列:"<<endl;for(i=0;i<n;i++){for(j=7;j<n;j++){if (fabs(qr[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<qr[i][j]<<"";}cout<<endl;}cout<<endl<<"矩阵A的全部特征值为"<<endl;for(i=0;i<n;i++){if(L[i][1]==0)cout<<"λ"<<i+1<<"="<<L[i][0]<<endl;elsecout<<"λ"<<i+1<<"="<<L[i][0]<<"+"<<L[i][1]<<"i"<<endl;}cout<<endl<<"实特征值对应的特征向量分别是："<<endl;for(k=0;k<n;k++){if(L[k][1]==0)//判断特征值是否为实特征值{for(i=0;i<n;i++)//重新录入矩阵A{for(j=0;j<n;j++){if(i==j)a[i][j]=1.5*cos(i+1+1.2*(j+1))-L[k][0];elsea[i][j]=sin(0.5*(i+1)+0.2*(j+1));}}Gauss(a,x);cout<<"λ"<<k+1<<"="<<L[k][0]<<""<<"对应的特征向量是:"<<endl;for(j=0;j<n;j++){cout<<x[j]<<endl;}}else continue;}}void Householder(double a[n][n])//拟上三角化函数{int i,j,k;int m=0;double d,c,h,t;double u[n],p[n],q[n],w[n];for(i=0;i<n-2;i++){for(j=i+2;j<n;j++){if(a[j][i]<=E)m=m+1;}if(m==(n-2-i))continue;for(j=i+1,d=0;j<n;j++){d=d+a[j][i]*a[j][i];}d=sqrt(d);c=-1*sgn(a[i+1][i])*d;h=c*c-c*a[i+1][i];for(j=i+2;j<n;j++){u[j]=a[j][i];}for(j=0;j<i+2;j++){u[j]=0;}u[i+1]=a[i+1][i]-c;for(j=0;j<n;j++){for(k=i+1,p[j]=0;k<n;k++){p[j]=a[k][j]*u[k]+p[j];}p[j]=p[j]/h;}for(j=0;j<n;j++){for(k=i+1,q[j]=0;k<n;k++){q[j]=a[j][k]*u[k]+q[j];}q[j]=q[j]/h;}for(j=0,t=0;j<n;j++){t=t+p[j]*u[j];}t=t/h;for(j=0;j<n;j++){w[j]=q[j]-t*u[j];}for(j=0;j<n;j++){for(k=0;k<n;k++){a[j][k]=a[j][k]-w[j]*u[k]-u[j]*p[k];}}}} //以上是拟上三角化的结果double sgn(double a)//符号函数{if(a>0)return 1;else if(a<0)return -1;return 0;}//以上是符号函数void QRfenjie(double a[n][n],double Q[n][n],double R[n][n])//矩阵的QR分解{int i,j,k;double d,c,h;double u[n],w[n],p[n];for(i=0;i<n;i++){for(j=0;j<n;j++){if (i==j) Q[i][j]=1; else Q[i][j]=0;}}for(i=0;i<n;i++){for(j=0;j<n;j++) R[i][j]=a[i][j];}for(i=0;i<n-1;i++){for(j=i,d=0;j<n;j++)d=d+R[j][i]*R[j][i];d=sqrt(d);c=-1*sgn(R[i][i])*d;h=c*c-c*R[i][i];for(j=i+1;j<n;j++) u[j]=R[j][i];for(j=0;j<i;j++) u[j]=0;u[i]=R[i][i]-c;for(j=0;j<n;j++)for(k=0,w[j]=0;k<n;k++) w[j]=Q[j][k]*u[k]+w[j];}for(j=0;j<n;j++){for(k=0;k<n;k++) Q[j][k]=Q[j][k]-w[j]*u[k]/h;}for(j=0;j<n;j++){for(k=i,p[j]=0;k<n;k++)p[j]=R[k][j]*u[k]+p[j];p[j]=p[j]/h;}for(j=0;j<n;j++){for(k=0;k<n;k++)R[j][k]=R[j][k]-u[j]*p[k];}}}//以上是矩阵的QR分解void MxM(double M[n][n],double A[n][n],double B[n][n],int m)//矩阵相乘{double C[n][n],D[n][n];int i,j,g;for(i=0;i<m;i++)for(j=0;j<m;j++){C[i][j]=A[i][j];D[i][j]=B[i][j];M[i][j]=0;}for(i=0;i<m;i++)for(j=0;j<m;j++)for(g=0;g<m;g++)M[i][j]=M[i][j]+C[i][g]*D[g][j];}void QR(double a[n][n],double L[n][2])//带双步位移的QR方法{double M[n][n],s,t,s1[2],s2[2],u[n],v[n],p[n],q[n],w[n],sum,h,c,d;int i,j,m,r,k;for(i=0;i<n;i++){for(j=0;j<2;j++){L[i][j]=0;}}m=n;for(k=1;k<=100;k++){if(fabs(a[m-1][m-2])<=E){L[m-1][0]=a[m-1][m-1];m=m-1;if(m==1){L[0][0]=a[0][0];break;}else if(m==0)break;elsecontinue;}else{solve(a,s1,s2,m);//调用解方程函数，将求解结果存到s1和s2中if(m==2){L[0][0]=s1[0];L[0][1]=s1[1];L[1][0]=s2[0];L[1][1]=s2[1];break;}else{if(fabs(a[m-2][m-3])<=E){L[m-2][0]=s1[0];L[m-2][1]=s1[1];L[m-1][0]=s2[0];L[m-1][1]=s2[1];m=m-2;if(m==1){L[0][0]=a[0][0];break;}else if(m==0)break;elsecontinue;}else{if(k==100){cout<<"未能得到全部特征值"<<endl;break;}else{t=a[m-1][m-1]*a[m-2][m-2]-a[m-1][m-2]*a[m-2][m-1];s=a[m-1][m-1]+a[m-2][m-2];MxM(M,a,a,m); //调用矩阵相乘函数构造矩阵Mfor(i=0;i<m;i++)for(j=0;j<m;j++)M[i][j]=M[i][j]-s*a[i][j];for(i=0;i<m;i++)M[i][i]=M[i][i]+t;for(r=0;r<(m-1);r++){for(i=(r+1);i<m;i++){if(M[i][r]==0)continue;else{sum=0;for(i=r;i<m;i++){sum+=pow(M[i][r],2);}d=sqrt(sum);if(M[r][r]!=0)c=-sgn(M[r][r])*d;else c=d;h=c*(c-M[r][r]);for(i=0;i<m;i++){if(i<r)u[i]=0;if(i==r)u[i]=M[r][r]-c;if(i>r)u[i]=M[i][r];}for(i=0;i<m;i++){sum=0;for(j=r;j<m;j++){sum+=M[j][i]*u[j];}v[i]=sum/h;}for(i=0;i<m;i++){for(j=0;j<m;j++){M[i][j]=M[i][j]-u[i]*v[j];} }for(i=0;i<m;i++){sum=0;for(j=r;j<m;j++){sum+=a[j][i]*u[j];}p[i]=sum/h;sum=0;for(j=r;j<m;j++){sum+=a[i][j]*u[j];}q[i]=sum/h;}sum=0;for(j=r;j<m;j++){sum+=p[j]*u[j];}t=sum/h;for(j=0;j<m;j++){w[j]=q[j]-t*u[j];}for(i=0;i<m;i++){for(j=0;j<m;j++){a[i][j]=a[i][j]-w[i]*u[j]-u[i]*p[j];}}}}}}}}}}for(i=0;i<n;i++){for(j=0;j<n;j++){if(a[i][j]<=E)a[i][j]=0;}}}void solve(double a[n][n],double s1[2],double s2[2],int m)//解方程函数{double s,t,det;t=a[m-1][m-1]*a[m-2][m-2]-a[m-1][m-2]*a[m-2][m-1];s=a[m-1][m-1]+a[m-2][m-2];det=s*s-4*t;if(det<0){det=-det;det=sqrt(det);s1[0]=s/2;s1[1]=det/2;s2[0]=s/2;s2[1]=-det/2;}else{det=sqrt(det);s1[0]=(s+det)/2;s1[1]=0;s2[0]=(s-det)/2;s2[1]=0;}}void Gauss(double a[n][n], double x[n])//定义列主元高斯消去法函数{int i,j,k,t;double sum,m[n][n],max,p;for(k=0;k<(n-1);k++){max=fabs(a[k][k]);t=k;for(i=(k+1);i<n;i++){if(max<fabs(a[i][k])){max=fabs(a[i][k]);t=i;}}if(t!=k){for(j=k;j<n;j++){p=a[k][j];a[k][j]=a[t][j];a[t][j]=p;}}for(i=(k+1);i<n;i++){m[i][k]=a[i][k]/a[k][k];for(j=(k+1);j<n;j++){a[i][j]=a[i][j]-m[i][k]*a[k][j];}}}x[n-1]=1;for(k=(n-2);k>=0;k--){sum=0;for(j=(k+1);j<n;j++){sum+=a[k][j]*x[j];}x[k]=(-sum)/a[k][k];}sum=0;for(i=0;i<n;i++){sum+=x[i]*x[i];}sum=sqrt(sum);for(i=0;i<n;i++){x[i]=x[i]/sum;}}三、处理结果四、结果讨论（1）此次编程涉及的函数众多，需要特别注意各函数之间的协调；（2）采用全局变量n节省了输入次数，简化编程；Q R矩阵，还需要另外的编辑函数来求解除Q和R，然后再求出Q R。