欢迎来到天天文库
浏览记录
ID:9982181
大小:34.00 KB
页数:3页
时间:2018-05-19
《雅克比过关法求特征值和特征向量》由会员上传分享,免费在线阅读,更多相关内容在行业资料-天天文库。
1、1.//////////////////////////////////////////////////////////////////////2.//求实对称矩阵特征值与特征向量的雅可比法3.//4.//参数:5.//1.doubledblEigenValue[]-一维数组,长度为矩阵的阶数,返回时存放特征值6.//2.CMatrix&mtxEigenVector-返回时存放特征向量矩阵,其中第i列为与7.//数组dblEigenValue中第j个特征值对应的特征向量8.//3.intnMaxIt-迭代次数,默认值为609.//4.doubleeps-计算精度,默
2、认值为0.00000110.//11.//返回值:BOOL型,求解是否成功12.//////////////////////////////////////////////////////////////////////13.BOOLCMatrix::JacobiEigenv(doubledblEigenValue[],CMatrix&mtxEigenVector,intnMaxIt/*=60*/,doubleeps/*=0.000001*/)14.{15.inti,j,p,q,u,w,t,s,l;16.doublefm,cn,sn,omega,x,y,d;17.1
3、8.if(!mtxEigenVector.Init(m_nNumColumns,m_nNumColumns))19.returnFALSE;20.21.l=1;22.for(i=0;i<=m_nNumColumns-1;i++)23.{24.mtxEigenVector.m_pData[i*m_nNumColumns+i]=1.0;25.for(j=0;j<=m_nNumColumns-1;j++)26.if(i!=j)27.mtxEigenVector.m_pData[i*m_nNumColumns+j]=0.0;//单位矩阵28.}29.30.while(TRU
4、E)31.{32.fm=0.0;33.for(i=1;i<=m_nNumColumns-1;i++)34.{35.for(j=0;j<=i-1;j++)36.{37.d=fabs(m_pData[i*m_nNumColumns+j]);38.if((i!=j)&&(d>fm))39.{40.fm=d;41.p=i;42.q=j;}//取绝对值最大的非对角线元素,并记住位置1.}2.}3.4.if(fm5、.}10.11.if(l>nMaxIt)12.returnFALSE;13.14.l=l+1;15.u=p*m_nNumColumns+q;16.w=p*m_nNumColumns+p;17.t=q*m_nNumColumns+p;18.s=q*m_nNumColumns+q;19.x=-m_pData;20.y=(m_pData[s]-m_pData[w])/2.0;21.omega=x/sqrt(x*x+y*y);22.23.if(y<0.0)24.omega=-omega;25.26.sn=1.0+sqrt(1.0-omega*omega);27.sn=ome6、ga/sqrt(2.0*sn);//sin28.cn=sqrt(1.0-sn*sn);//cos29.fm=m_pData[w];30.m_pData[w]=fm*cn*cn+m_pData[s]*sn*sn+m_pData*omega;31.m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData*omega;32.m_pData=0.0;33.m_pData[t]=0.0;34.for(j=0;j<=m_nNumColumns-1;j++)35.{36.if((j!=p)&&(j!=q))37.{38.u=p*m_nNumColu7、mns+j;w=q*m_nNumColumns+j;39.fm=m_pData;40.m_pData=fm*cn+m_pData[w]*sn;41.m_pData[w]=-fm*sn+m_pData[w]*cn;42.}43.}1.for(i=0;i<=m_nNumColumns-1;i++)2.{3.if((i!=p)&&(i!=q))4.{5.u=i*m_nNumColumns+p;6.w=i*m_nNumColumns+q;7.fm=m_pData;8.m_pData=fm*cn+m_pData[w]*sn;9.m_pData[w]=-fm*sn+m_pD
5、.}10.11.if(l>nMaxIt)12.returnFALSE;13.14.l=l+1;15.u=p*m_nNumColumns+q;16.w=p*m_nNumColumns+p;17.t=q*m_nNumColumns+p;18.s=q*m_nNumColumns+q;19.x=-m_pData;20.y=(m_pData[s]-m_pData[w])/2.0;21.omega=x/sqrt(x*x+y*y);22.23.if(y<0.0)24.omega=-omega;25.26.sn=1.0+sqrt(1.0-omega*omega);27.sn=ome
6、ga/sqrt(2.0*sn);//sin28.cn=sqrt(1.0-sn*sn);//cos29.fm=m_pData[w];30.m_pData[w]=fm*cn*cn+m_pData[s]*sn*sn+m_pData*omega;31.m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData*omega;32.m_pData=0.0;33.m_pData[t]=0.0;34.for(j=0;j<=m_nNumColumns-1;j++)35.{36.if((j!=p)&&(j!=q))37.{38.u=p*m_nNumColu
7、mns+j;w=q*m_nNumColumns+j;39.fm=m_pData;40.m_pData=fm*cn+m_pData[w]*sn;41.m_pData[w]=-fm*sn+m_pData[w]*cn;42.}43.}1.for(i=0;i<=m_nNumColumns-1;i++)2.{3.if((i!=p)&&(i!=q))4.{5.u=i*m_nNumColumns+p;6.w=i*m_nNumColumns+q;7.fm=m_pData;8.m_pData=fm*cn+m_pData[w]*sn;9.m_pData[w]=-fm*sn+m_pD
此文档下载收益归作者所有