特征向量求法2.docx

/数值计算程序-特征值和特征向量 / /约化对称矩阵为三对角对称矩阵 /利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵 /a-长度为n*n的数组，存放n阶实对称矩阵 /n-矩阵的阶数 /q-长度为n*n的数组，返回时存放Householder变换矩阵 /b-长度为n的数组，返回时存放三对角阵的主对角线元素 /c-长度为n的数组，返回时前n-1个元素存放次对角线元素 void eastrq(double a,int n,double q,double b,double c); / /求实对称三对角对称矩阵的全部特征值及特征向量 /利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量 /n-矩阵的阶数 /b-长度为n的数组，返回时存放三对角阵的主对角线元素 /c-长度为n的数组，返回时前n-1个元素存放次对角线元素 /q-长度为n*n的数组，若存放单位矩阵，则返回实对称三对角矩阵的特征向量组 / 若存放Householder变换矩阵，则返回实对称矩阵A的特征向量组 /a-长度为n*n的数组，存放n阶实对称矩阵 int ebstq(int n,double b,double c,double q,double eps,int l); / /约化实矩阵为赫申伯格(Hessen berg)矩阵 /利用初等相似变换将n阶实矩阵约化为上H矩阵 /a-长度为n*n的数组，存放n阶实矩阵,返回时存放上H矩阵 /n-矩阵的阶数 void echbg(double a,int n); / /求赫申伯格(Hessen berg)矩阵的全部特征值 /返回值小于0表示超过迭代jt次仍未达到精度要求 /返回值大于0表示正常返回 /利用带原点位移的双重步QR方法求上H矩阵的全部特征值 /a-长度为n*n的数组，存放上H矩阵 /n-矩阵的阶数 /u-长度为n的数组，返回n个特征值的实部 /v-长度为n的数组，返回n个特征值的虚部 /eps-控制精度要求 /jt-整型变量，控制最大迭代次数 int edqr(double a,int n,double u,double v,double eps,int jt); / /求实对称矩阵的特征值及特征向量的雅格比法 /利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量 /返回值小于0表示超过迭代jt次仍未达到精度要求 /返回值大于0表示正常返回 /a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值 /n-矩阵的阶数 /u-长度为n*n的数组，返回特征向量(按列存储) /eps-控制精度要求 /jt-整型变量，控制最大迭代次数 int eejcb(double a,int n,double v,double eps,int jt); / 选自 每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。 今天都给贴出来了 #include stdio.h #include math.h /约化对称矩阵为三对角对称矩阵 /利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵 /a-长度为n*n的数组，存放n阶实对称矩阵 /n-矩阵的阶数 /q-长度为n*n的数组，返回时存放Householder变换矩阵 /b-长度为n的数组，返回时存放三对角阵的主对角线元素 /c-长度为n的数组，返回时前n-1个元素存放次对角线元素 void eastrq(double a,int n,double q,double b,double c) int i,j,k,u,v; double h,f,g,h2; for (i=0; i=n-1; i+) for (j=0; j=1; i-) h=0.0; if (i1) for (k=0; k0.0) ci-1=-ci-1; h=h-qu*ci-1; qu=qu-ci-1; f=0.0; for (j=0; j=i-1; j+) qj*n+i=qi*n+j/h; g=0.0; for (k=0; k=j; k+) g=g+qj*n+k*qi*n+k; if (j+1=i-1) for (k=j+1; k=i-1; k+) g=g+qk*n+j*qi*n+k; cj-1=g/h; f=f+g*qj*n+i; h2=f/(h+h); for (j=0; j=i-1; j+) f=qi*n+j; g=cj-1-h2*f; cj-1=g; for (k=0; k=j; k+) u=j*n+k; qu=qu-f*ck-1-g*qi*n+k; bi=h; b0=0.0; for (i=0; i=0) for (j=0; j=i-1; j+) g=0.0; for (k=0; k=i-1; k+) g=g+qi*n+k*qk*n+j; for (k=0; k=0) for (j=0; j=i-1; j+) qi*n+j=0.0; qj*n+i=0.0; return; /求实对称三对角对称矩阵的全部特征值及特征向量 /利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量 /n-矩阵的阶数 /b-长度为n的数组，返回时存放三对角阵的主对角线元素 /c-长度为n的数组，返回时前n-1个元素存放次对角线元素 /q-长度为n*n的数组，若存放单位矩阵，则返回实对称三对角矩阵的特征向量组 / 若存放Householder变换矩阵，则返回实对称矩阵A的特征向量组 /a-长度为n*n的数组，存放n阶实对称矩阵 int ebstq(int n,double b,double c,double q,double eps,int l) int i,j,k,m,it,u,v; double d,f,h,g,p,r,e,s; cn-1=0.0; d=0.0; f=0.0; for (j=0; jd) d=h; m=j; while (md) m=m+1; if (m!=j) do if (it=l) printf(failn); return(-1); it=it+1; g=bj; p=(bj+1-g)/(2.0*cj); r=sqrt(p*p+1.0); if (p=0.0) bj=cj/(p+r); else bj=cj/(p-r); h=g-bj; for (i=j+1; i=j; i-) g=e*ci; h=e*p; if (fabs(p)=fabs(ci) e=ci/p; r=sqrt(e*e+1.0); ci+1=s*p*r; s=e/r; e=1.0/r; else e=p/ci; r=sqrt(e*e+1.0); ci+1=s*ci*r; s=1.0/r; e=e/r; p=e*bi-s*g; bi+1=h+s*(e*g+s*bi); for (k=0; kd); bj=bj+f; for (i=0; i=n-1; i+) k=i; p=bi; if (i+1=n-1) j=i+1; while (j=n-1)&(bj=p) k=j; p=bj; j=j+1; if (k!=i) bk=bi; bi=p; for (j=0; j=n-1; j+) u=j*n+i; v=j*n+k; p=qu; qu=qv; qv=p; return(1); /约化实矩阵为赫申伯格(Hessen berg)矩阵 /利用初等相似变换将n阶实矩阵约化为上H矩阵 /a-长度为n*n的数组，存放n阶实矩阵,返回时存放上H矩阵 /n-矩阵的阶数 void echbg(double a,int n) int i,j,k,u,v; double d,t; for (k=1; k=n-2; k+) d=0.0; for (j=k; jfabs(d) d=t; i=j; if (fabs(d)+1.0!=1.0) if (i!=k) for (j=k-1; j=n-1; j+) u=i*n+j; v=k*n+j; t=au; au=av; av=t; for (j=0; j=n-1; j+) u=j*n+i; v=j*n+k; t=au; au=av; av=t; for (i=k+1; i=n-1; i+) u=i*n+k-1; t=au/d; au=0.0; for (j=k; j=n-1; j+) v=i*n+j; av=av-t*ak*n+j; for (j=0; j0)&(fabs(al*n+l-1)eps*(fabs(a(l-1)*n+l-1)+fabs(al*n+l) l=l-1; ii=(m-1)*n+m-1; jj=(m-1)*n+m-2; kk=(m-2)*n+m-1; ll=(m-2)*n+m-2; if (l=m-1) um-1=a(m-1)*n+m-1; vm-1=0.0; m=m-1; it=0; else if (l=m-2) b=-(aii+all); c=aii*all-ajj*akk; w=b*b-4.0*c; y=sqrt(fabs(w); if (w0.0) xy=1.0; if (b=jt) printf(failn); return(-1); it=it+1; for (j=l+2; j=m-1; j+) aj*n+j-2=0.0; for (j=l+3; j=m-1; j+) aj*n+j-3=0.0; for (k=l; k=m-2; k+) if (k!=l) p=ak*n+k-1; q=a(k+1)*n+k-1; r=0.0; if (k!=m-2) r=a(k+2)*n+k-1; else x=aii+all; y=all*aii-akk*ajj; ii=l*n+l; jj=l*n+l+1; kk=(l+1)*n+l; ll=(l+1)*n+l+1; p=aii*(aii-x)+ajj*akk+y; q=akk*(aii+all-x); r=akk*a(l+2)*n+l+1; if (fabs(p)+fabs(q)+fabs(r)!=0.0) xy=1.0; if (p0.0) xy=-1.0; s=xy*sqrt(p*p+q*q+r*r); if (k!=l) ak*n+k-1=-s; e=-q/s; f=-r/s; x=-p/s; y=-x-f*r/(p+s); g=e*r/(p+s); z=-x-e*q/(p+s); for (j=k; j=m-1) j=m-1; for (i=l; i=j; i+) ii=i*n+k; jj=i*n+k+1; p=x*aii+e*ajj; q=e*aii+y*ajj; r=f*aii+g*ajj; if (k!=m-2) kk=i*n+k+2; p=p+f*akk; q=q+g*akk; r=r+z*akk; akk=r; ajj=q; aii=p; return(1); /求实对称矩阵的特征值及特征向量的雅格比法 /利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量 /返回值小于0表示超过迭代jt次仍未达到精度要求 /返回值大于0表示正常返回 /a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值 /n-矩阵的阶数 /u-长度为n*n的数组，返回特征向量(按列存储) /eps-控制精度要求 /jt-整型变量，控制最大迭代次数 int eejcb(double a,int n,double v,double eps,int jt) int i,j,p,q,u,w,t,s,l; double fm,cn,sn,omega,x,y,d; l=1; for (i=0; i=n-1; i+) vi*n+i=1.0; for (j=0; j=n-1; j+) if (i!=j) vi*n+j=0.0; while (1=1) fm=0.0; for (i=0; i=n-1; i+) for (j=0; jfm) fm=d; p=i; q=j; if (fmjt) return(-1); l=l+1; u=p*n+q; w=p*n+p; t=q*n+p; s=q*n+q; x=-au; y=(as-aw)/2.0; omega=x/sqrt(x*x+y*y); if (y0.0) omega=-omega; sn=1.0+sqrt(1.0-omega*omega); sn=omega/sqrt(2.0*sn); cn=sqrt(1.0-sn*sn); fm=aw; aw=fm*cn*cn+as*sn*sn+au*omega; as=fm*sn*sn+as*cn*cn-au*omega; au=0.0; at=0.0; for (j=