数值分析计算实习二

1、数值分析计算实习二一 算法描述：用QR方法求矩阵的特征值和特征向量1.矩阵的拟上三角化实矩阵可以通过构造household矩阵的方法相似变换化为拟上三角矩阵。为节省计算量，使用教材上的算法。对完成拟上三角化的矩阵进行QR分解能减少计算量。2.带双步位移的QR方法QR方法产生的矩阵序列本质上收敛于分块上三角阵。对角块均是一阶或二阶子块，对每个子块求特征值就得到整个矩阵的所有特征值。为加速收敛，采用带双步位移的QR分解法。记最大迭代次数为Max，取为1000，迭代精度为e=1e-12。每次迭代前检查，若下一个子块满足条件（一阶或二阶），输所有特征值，计算完成。若达到最大迭代次数仍未输出所有特征值，

2、则返回错误信息，计算失败。采用书本上的算法可以节省计算量。3.用高斯消元法解特征向量。对每个实特征值i，解十阶方程组A-iIX=0,因为系数矩阵非满秩的，统一置特征向量的最后一个值为1。二 C语言源代码：/带双步位移的QR方法#include#include#include#define e 1e-12 /迭代的精度水平#define Max 1000 /迭代的最大次数typedef structdouble Re;double Im;complex;void creatematrix(double *a); /创建题中要求的矩阵int hessenberg(double *a); /对矩阵进

3、行拟上三角化int muiltiply(double *a,double *b,double *c,int m); /计算m阶矩阵的乘法a*b=cint printmatrix(double *a); /打印矩阵void QRmethod(double *a); /带双步位移的QR分解法void check(double *a); /将满足精度要求可看作零的元素全部置为0void output(complex *L); /输出特征值和特征向量void calculate(double *M, double *a,int m); /执行QR方法中迭代计算的步骤void gauss(double

4、lambda); /高斯消元void maxline(double *ptr, int n, int k);int main()double *a;int i;a=(double *)malloc(10*sizeof(double *);for(i=0;i10;i+)ai=(double *)malloc(10*sizeof(double);creatematrix(a);hessenberg(a);check(a);printmatrix(a);/打印已完成上三角化的矩阵QRmethod(a);void creatematrix(double *a)int i,j;for(i=0;i10;i

5、+)for(j=0;j10;j+)if (i!=j)aij=sin(0.5*(i+1)+0.2*(j+1);elseaij=1.52*cos(i+1+1.2*(j+1);int hessenberg(double *a) /构造的拟上三角阵存储在原来a的地址中int r,i,j;double c,d,h,t,u10,p10,q10,w10;for(r=0;r8;r+)check(a);c=0;d=0;h=0;for(i=r+2;i10;i+)d+=air*air;if(d=0)continue;elsed+=ar+1r*ar+1r;d=pow(d,0.5);if(ar+1r!=0)c=-fab

6、s(ar+1r)/ar+1r*d;else c=d;h=c*c-c*ar+1r;for(i=0;ir+1;i+)ui=0;ur+1=ar+1r-c;for(i=r+2;i10;i+) ui=air;for(i=0;i10;i+)pi=0;qi=0;for(j=0;j10;j+)pi+=aji*uj/h;qi+=aij*uj/h;t=0;for(i=0;i10;i+)t+=pi*ui/h;for(i=0;i10;i+)wi=qi-t*ui;for(i=0;i10;i+)for(j=0;j10;j+)aij-=(wi*uj+ui*pj);int muiltiply(double *a,double

7、 *b,double *c,int m)double *ctemp;/临时存储计算结果；int i,j,k;ctemp=(double *)malloc(m*sizeof(double *);for(i=0;im;i+)ctempi=(double*)malloc(m*sizeof(double);for(i=0;im;i+)for(j=0;jm;j+)ctempij=0;for(k=0;km;k+)ctempij+=aik*bkj;for(i=0;im;i+)for(j=0;jm;j+)cij=ctempij;for(i=0;im;i+)free(ctempi);free(ctemp);re

8、turn 0;int printmatrix(double *a)int i,j;for(i=0;i10;i+)for(j=0;j10;j+)printf(%.12e ,aij);printf(n);printf(nnn);return 0;void QRmethod(double *a)int k,m,i,j,r;double s,t,det;double *M;complex L10;/为Q，R，M分配地址M=(double*)malloc(10*sizeof(double *);for(k=0;k10;k+)Mk=(double*)malloc(10*sizeof(double);m=9

9、;r=0;for(k=0;kMax;k+)if(m=0)Lr.Re=amm; Lr.Im=0;break;else if(m0) break;if(fabs(amm-1)0)Lr.Re=(amm+am-1m-1)/2+sqrt(det)/2; Lr.Im=0; Lr+1.Re=(amm+am-1m-1)/2-sqrt(det)/2; Lr+1.Im=0;else Lr.Re=(amm+am-1m-1)/2; Lr.Im=sqrt(-det)/2; Lr+1.Re=(amm+am-1m-1)/2;Lr+1.Im=-sqrt(-det)/2;m-=2;r+=2;continue;else if(f

10、abs(am-1m-2)0)Lr.Re=(amm+am-1m-1)/2+sqrt(det)/2; Lr.Im=0; Lr+1.Re=(amm+am-1m-1)/2-sqrt(det)/2; Lr+1.Im=0;else Lr.Re=(amm+am-1m-1)/2; Lr.Im=sqrt(-det)/2; Lr+1.Re=(amm+am-1m-1)/2;Lr+1.Im=-sqrt(-det)/2;m-=2;r+=2;continue;elses=am-1m-1+amm;t=am-1m-1*amm-amm-1*am-1m;muiltiply(a,a,M,m+1);for(i=0;i10;i+)fo

11、r(j=0;j10;j+)Mij-=s*aij;Mii+=t;calculate(M,a,m+1);check(a);check(a);printmatrix(a);output(L);void check(double *a)int i,j;for(i=0;i10;i+)for(j=0;j10;j+)if(aij-e)aij=0;void output(complex *L)int r;for(r=0;r10;r+)if(Lr.Im=0)printf(lambda%d=%en,r+1,Lr.Re);gauss(Lr.Re);else printf(lambda%d=%e+i*%en,r+1,

12、Lr.Re,Lr.Im);void calculate(double *M,double *a,int m)int r,i,j;double c,d,h,t,u10,v10,p10,q10,w10;for(r=0;rm-1;r+)check(M);c=0;d=0;h=0;for(i=r+1;im;i+)d+=Mir*Mir;if(fabs(d)=0)continue;elsed+=Mrr*Mrr;d=pow(d,0.5);if(Mrr!=0)c=-fabs(Mrr)/Mrr*d;else c=d;h=c*c-c*Mrr;for(i=0;ir;i+)ui=0;ur=Mrr-c;for(i=r+1

13、;im;i+) ui=Mir;for(i=0;im;i+)vi=0;for(j=0;jm;j+)vi+=Mji*uj/h;for(i=0;im;i+)pi=0;qi=0;for(j=0;jm;j+)Mij-=ui*vj;pi+=aji*uj/h;qi+=aij*uj/h;t=0;for(i=0;im;i+)t+=pi*ui/h;for(i=0;im;i+)wi=qi-t*ui;for(i=0;im;i+)for(j=0;jm;j+)aij-=(wi*uj+ui*pj);void gauss(double lambda)double *a;double *X;double m,sigma;int

14、 i,j,k,n=10;X=(double *)malloc(10*sizeof(double);a=(double *)malloc(10*sizeof(double *);for(i=0;in;i+)ai=(double *)malloc(11)*sizeof(double);creatematrix(a);for(i=0;in;i+)aii-=lambda; ai10=0;for(k=0;kn-1;k+)maxline(a,n,k);for(i=k+1;in;i+)m=aik/akk;for(j=k;j=0;k-)sigma=0;for(j=k+1;jn;j+)sigma+=akj*Xj

15、;Xk=(akn-sigma)/akk;printf(对应的特征向量：( );for(i=0;in;i+)printf(%lf ,Xi);printf()n);for(i=0;i10;i+)free(ai);free(a);free(X);void maxline(double *ptr, int n, int k)double c;int i,M;M=k;for(i=k;ifabs(ptrMk)M=i;if(Mk) for(i=k;in+1;i+) c=ptrki; ptrki=ptrMi; ptrMi=c;四计算结果：An-1=-8.1e-001 -9.6e-002 -1.7e+000 -

16、7.7e-001 1.6e-001 -1.9e+000 -8.6e-002 9.3e-001 -6.8e-001 1.5e-001-2.6e+000 2.8e+000 1.6e+000 3.4e-001 2.5e-001 2.1e+000 1.9e-001 -1.0e+000 6.9e-001 -4.0e-0010.0e+000 1.6e+000 -1.7e+000 -1.8e+000 -6.4e-001 -1.2e+000 2.0e-001 6.5e-001 -1.2e-001 3.5e-0010.0e+000 0.0e+000 -1.4e+000 -1.2e+000 1.4e+000 -1

17、.5e+000 1.6e-001 4.1e-001 -1.5e-001 1.5e-0010.0e+000 0.0e+000 0.0e+000 1.2e+0008.0e-001 4.8e-001 -4.9e-002 -4.8e-001 2.4e-001 -1.5e-0010.0e+000 0.0e+000 0.0e+000 0.0e+000-7.9e-001 -1.6e+000 -2.7e-001 -2.1e-001 6.0e-001 2.7e-0010.0e+000 0.0e+000 0.0e+000 0.0e+0000.0e+000 -7.2e-001 -7.9e-003 9.7e-001

18、-1.4e-001 2.1e-0020.0e+000 0.0e+000 0.0e+000 0.0e+0000.0e+000 0.0e+000 7.6e-001 -4.6e-001 5.5e-001 -1.2e-0010.0e+000 0.0e+000 0.0e+000 0.0e+0000.0e+000 0.0e+000 0.0e+000 7.4e-0011.8e-001 -3.3e-0010.0e+000 0.0e+000 0.0e+000 0.0e+0000.0e+000 0.0e+000 0.0e+000 0.0e+000-4.4e-001 4.9e-001QR迭代完成的矩阵：3.6e+0

19、00 8.9e-001 -9.1e-001 -8.9e-002 2.4e-001 9.4e-001 -1.2e+000 -7.9e-002 9.2e-001 -4.6e-0010.0e+000 -2.6e+000 -2.6e+000 4.1e-002 -4.5e-002 1.1e+000 -1.0e+000 1.9e-001 5.0e-001 6.0e-0020.0e+000 3.5e-001 -2.1e+000 6.6e-001 2.1e-002 1.2e+000 -7.0e-001 3.9e-001 -1.9e+000 1.9e+0000.0e+000 0.0e+000 0.0e+000

20、1.6e+0001.6e-002 -5.0e-001 4.1e-001 -7.8e-003 8.8e-002 -1.2e-0010.0e+000 0.0e+000 0.0e+000 -1.4e-006 -1.5e+000 -7.9e-002 8.6e-002 -3.4e-002 -1.0e-001 -1.9e-0010.0e+000 0.0e+000 0.0e+000 0.0e+0000.0e+000 -9.6e-001 -6.6e-001 -2.3e-001 5.2e-001 -5.8e-0010.0e+000 0.0e+000 0.0e+000 0.0e+0000.0e+000 1.9e-002 -9.1e-001 -1.2e-001 -2.3e-001 -4.6e-0020.0e+000 0.0e+000 0.0e+000 0.0e+0000.0e+000 0.0e+000 0.0e+000 9.9e-0011.3e-001 -1.0e-0010.0e+000 0.0e+000 0.0e+000 0.0e+0000.0e+000 0.0e+000 0.0e+000 0.0e+0006.1e-001 -2.2e-0010.0e+000 0.0e+000 0.0e+000 0.0