算法及编程语言 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

风花雪月 发表于 2005-9-20 10:07

[分享]求矩阵特征值和特征向量的一个小程序

代码较长，如果不能执行，就是要建立结构体，大家试试吧，希望能用。

//////////////////////////////////////////////////////////////////////
// 实对称三对角阵的全部特征值与特征向量的计算
//
// 参数：
// 1. double dblB[] - 一维数组，长度为矩阵的阶数，传入对称三对角阵的主对角线元素；
// 返回时存放全部特征值。
// 2. double dblC[] - 一维数组，长度为矩阵的阶数，前n-1个元素传入对称三对角阵的次对角线元素
// 3. CMatrix& mtxQ - 如果传入单位矩阵，则返回实对称三对角阵的特征值向量矩阵；
// 如果传入MakeSymTri函数求得的矩阵A的豪斯荷尔德变换的乘积矩阵Q，则返回矩阵A的
// 特征值向量矩阵。其中第i列为与数组dblB中第j个特征值对应的特征向量。
// 4. int nMaxIt - 迭代次数，默认值为60
// 5. double eps - 计算精度，默认值为0.000001
//
// 返回值：BOOL型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::SymTriEigenv(double dblB[], double dblC[], CMatrix& mtxQ, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{
int i,j,k,m,it,u,v;
double d,f,h,g,p,r,e,s;

// 初值
int n = mtxQ.GetNumColumns();
dblC=0.0;
d=0.0;
f=0.0;

// 迭代计算

for (j=0; j<=n-1; j++)
{
it=0;
h=eps*(fabs(dblB)+fabs(dblC));
if (h>d)
d=h;

m=j;
while ((m<=n-1)&&(fabs(dblC)>d))
m=m+1;

if (m!=j)
{
do
{
if (it==nMaxIt)
return FALSE;

it=it+1;
g=dblB;
p=(dblB-g)/(2.0*dblC);
r=sqrt(p*p+1.0);
if (p>=0.0)
dblB=dblC/(p+r);
else
dblB=dblC/(p-r);

h=g-dblB;
for (i=j+1; i<=n-1; i++)
dblB=dblB-h;

f=f+h;
p=dblB;
e=1.0;
s=0.0;
for (i=m-1; i>=j; i--)
{
g=e*dblC;
h=e*p;
if (fabs(p)>=fabs(dblC))
{
e=dblC/p;
r=sqrt(e*e+1.0);
dblC=s*p*r;
s=e/r;
e=1.0/r;
}
else
{
e=p/dblC;
r=sqrt(e*e+1.0);
dblC=s*dblC*r;
s=1.0/r;
e=e/r;
}

p=e*dblB-s*g;
dblB=h+s*(e*g+s*dblB);
for (k=0; k<=n-1; k++)
{
u=k*n+i+1;
v=u-1;
h=mtxQ.m_pData;
mtxQ.m_pData=s*mtxQ.m_pData+e*h;
mtxQ.m_pData=e*mtxQ.m_pData-s*h;
}
}

dblC=s*p;
dblB=e*p;

} while (fabs(dblC)>d);
}

dblB=dblB+f;
}

for (i=0; i<=n-1; i++)
{
k=i;
p=dblB;
if (i+1<=n-1)
{
j=i+1;
while ((j<=n-1)&&(dblB<=p))
{
k=j;
p=dblB;
j=j+1;
}
}

if (k!=i)
{
dblB=dblB;
dblB=p;
for (j=0; j<=n-1; j++)
{
u=j*n+i;
v=j*n+k;
p=mtxQ.m_pData;
mtxQ.m_pData=mtxQ.m_pData;
mtxQ.m_pData=p;
}
}
}

return TRUE;
}

//////////////////////////////////////////////////////////////////////
// 约化一般实矩阵为赫申伯格矩阵的初等相似变换法
//
// 参数：无
//
// 返回值：无
//////////////////////////////////////////////////////////////////////
void CMatrix::MakeHberg()
{
int i,j,k,u,v;
double d,t;

for (k=1; k<=m_nNumColumns-2; k++)
{
d=0.0;
for (j=k; j<=m_nNumColumns-1; j++)
{
u=j*m_nNumColumns+k-1;
t=m_pData;
if (fabs(t)>fabs(d))
{
d=t;
i=j;
}
}

if (d != 0.0)
{
if (i!=k)
{
for (j=k-1; j<=m_nNumColumns-1; j++)
{
u=i*m_nNumColumns+j;
v=k*m_nNumColumns+j;
t=m_pData;
m_pData=m_pData;
m_pData=t;
}

for (j=0; j<=m_nNumColumns-1; j++)
{
u=j*m_nNumColumns+i;
v=j*m_nNumColumns+k;
t=m_pData;
m_pData=m_pData;
m_pData=t;
}
}

for (i=k+1; i<=m_nNumColumns-1; i++)
{
u=i*m_nNumColumns+k-1;
t=m_pData/d;
m_pData=0.0;
for (j=k; j<=m_nNumColumns-1; j++)
{
v=i*m_nNumColumns+j;
m_pData=m_pData-t*m_pData;
}

for (j=0; j<=m_nNumColumns-1; j++)
{
v=j*m_nNumColumns+k;
m_pData=m_pData+t*m_pData;
}
}
}
}
}

[ 本帖最后由风花雪月于 2006-8-24 03:31 编辑 ]

风花雪月 发表于 2005-9-20 10:07

回复：（风花雪月）[分享]求矩阵特征值和特征向量的一...

//////////////////////////////////////////////////////////////////////
// 求赫申伯格矩阵全部特征值的QR方法
//
// 参数：
// 1. double dblU[] - 一维数组，长度为矩阵的阶数，返回时存放特征值的实部
// 2. double dblV[] - 一维数组，长度为矩阵的阶数，返回时存放特征值的虚部
// 3. int nMaxIt - 迭代次数，默认值为60
// 4. double eps - 计算精度，默认值为0.000001
//
// 返回值：BOOL型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::HBergEigenv(double dblU[], double dblV[], int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{
int m,it,i,j,k,l,ii,jj,kk,ll;
double b,c,w,g,xy,p,q,r,x,s,e,f,z,y;

int n = m_nNumColumns;

it=0;
m=n;
while (m!=0)
{
l=m-1;
while ((l>0)&&(fabs(m_pData) >
eps*(fabs(m_pData[(l-1)*n+l-1])+fabs(m_pData))))
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)
{
dblU=m_pData[(m-1)*n+m-1];
dblV=0.0;
m=m-1;
it=0;
}
else if (l==m-2)
{
b=-(m_pData+m_pData);
c=m_pData*m_pData-m_pData*m_pData;
w=b*b-4.0*c;
y=sqrt(fabs(w));
if (w>0.0)
{
xy=1.0;
if (b<0.0)
xy=-1.0;
dblU=(-b-xy*y)/2.0;
dblU=c/dblU;
dblV=0.0; dblV=0.0;
}
else
{
dblU=-b/2.0;
dblU=dblU;
dblV=y/2.0;
dblV=-dblV;
}

m=m-2;
it=0;
}
else
{
if (it>=nMaxIt)
return FALSE;

it=it+1;
for (j=l+2; j<=m-1; j++)
m_pData=0.0;
for (j=l+3; j<=m-1; j++)
m_pData=0.0;
for (k=l; k<=m-2; k++)
{
if (k!=l)
{
p=m_pData;
q=m_pData[(k+1)*n+k-1];
r=0.0;
if (k!=m-2)
r=m_pData[(k+2)*n+k-1];
}
else
{
x=m_pData+m_pData;
y=m_pData*m_pData-m_pData*m_pData;
ii=l*n+l;
jj=l*n+l+1;
kk=(l+1)*n+l;
ll=(l+1)*n+l+1;
p=m_pData*(m_pData-x)+m_pData*m_pData+y;
q=m_pData*(m_pData+m_pData-x);
r=m_pData*m_pData[(l+2)*n+l+1];
}

if ((fabs(p)+fabs(q)+fabs(r))!=0.0)
{
xy=1.0;
if (p<0.0)
xy=-1.0;
s=xy*sqrt(p*p+q*q+r*r);
if (k!=l)
m_pData=-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++)
{
ii=k*n+j;
jj=(k+1)*n+j;
p=x*m_pData+e*m_pData;
q=e*m_pData+y*m_pData;
r=f*m_pData+g*m_pData;
if (k!=m-2)
{
kk=(k+2)*n+j;
p=p+f*m_pData;
q=q+g*m_pData;
r=r+z*m_pData;
m_pData=r;
}

m_pData=q; m_pData=p;
}

j=k+3;
if (j>=m-1)
j=m-1;

for (i=l; i<=j; i++)
{
ii=i*n+k;
jj=i*n+k+1;
p=x*m_pData+e*m_pData;
q=e*m_pData+y*m_pData;
r=f*m_pData+g*m_pData;
if (k!=m-2)
{
kk=i*n+k+2;
p=p+f*m_pData;
q=q+g*m_pData;
r=r+z*m_pData;
m_pData=r;
}

m_pData=q;
m_pData=p;
}
}
}
}
}

return TRUE;
}

//////////////////////////////////////////////////////////////////////
// 求实对称矩阵特征值与特征向量的雅可比法
//
// 参数：
// 1. double dblEigenValue[] - 一维数组，长度为矩阵的阶数，返回时存放特征值
// 2. CMatrix& mtxEigenVector - 返回时存放特征向量矩阵，其中第i列为与
// 数组dblEigenValue中第j个特征值对应的特征向量
// 3. int nMaxIt - 迭代次数，默认值为60
// 4. double eps - 计算精度，默认值为0.000001
//
// 返回值：BOOL型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::JacobiEigenv(double dblEigenValue[], CMatrix& mtxEigenVector, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{
int i,j,p,q,u,w,t,s,l;
double fm,cn,sn,omega,x,y,d;

if (! mtxEigenVector.Init(m_nNumColumns, m_nNumColumns))
return FALSE;

l=1;
for (i=0; i<=m_nNumColumns-1; i++)
{
mtxEigenVector.m_pData=1.0;
for (j=0; j<=m_nNumColumns-1; j++)
if (i!=j)
mtxEigenVector.m_pData=0.0;
}

while (TRUE)
{
fm=0.0;
for (i=1; i<=m_nNumColumns-1; i++)
{
for (j=0; j<=i-1; j++)
{
d=fabs(m_pData);
if ((i!=j)&&(d>fm))
{
fm=d;
p=i;
q=j;
}
}
}

if (fm<eps)
{
for (i=0; i<m_nNumColumns; ++i)
dblEigenValue = GetElement(i,i);
return TRUE;
}

if (l>nMaxIt)
return FALSE;

l=l+1;
u=p*m_nNumColumns+q;
w=p*m_nNumColumns+p;
t=q*m_nNumColumns+p;
s=q*m_nNumColumns+q;
x=-m_pData;
y=(m_pData-m_pData)/2.0;
omega=x/sqrt(x*x+y*y);

if (y<0.0)
omega=-omega;

sn=1.0+sqrt(1.0-omega*omega);
sn=omega/sqrt(2.0*sn);
cn=sqrt(1.0-sn*sn);
fm=m_pData;
m_pData=fm*cn*cn+m_pData*sn*sn+m_pData*omega;
m_pData=fm*sn*sn+m_pData*cn*cn-m_pData*omega;
m_pData=0.0;
m_pData=0.0;
for (j=0; j<=m_nNumColumns-1; j++)
{
if ((j!=p)&&(j!=q))
{
u=p*m_nNumColumns+j; w=q*m_nNumColumns+j;
fm=m_pData;
m_pData=fm*cn+m_pData*sn;
m_pData=-fm*sn+m_pData*cn;
}
}

for (i=0; i<=m_nNumColumns-1; i++)
{
if ((i!=p)&&(i!=q))
{
u=i*m_nNumColumns+p;
w=i*m_nNumColumns+q;
fm=m_pData;
m_pData=fm*cn+m_pData*sn;
m_pData=-fm*sn+m_pData*cn;
}
}

for (i=0; i<=m_nNumColumns-1; i++)
{
u=i*m_nNumColumns+p;
w=i*m_nNumColumns+q;
fm=mtxEigenVector.m_pData;
mtxEigenVector.m_pData=fm*cn+mtxEigenVector.m_pData*sn;
mtxEigenVector.m_pData=-fm*sn+mtxEigenVector.m_pData*cn;
}
}

for (i=0; i<m_nNumColumns; ++i)
dblEigenValue = GetElement(i,i);

return TRUE;
}

[ 本帖最后由风花雪月于 2006-11-14 19:23 编辑 ]

风花雪月 发表于 2005-9-20 10:07

回复：（风花雪月）[分享]求矩阵特征值和特征向量的一...

//////////////////////////////////////////////////////////////////////
// 求实对称矩阵特征值与特征向量的雅可比过关法
//
// 参数：
// 1. double dblEigenValue[] - 一维数组，长度为矩阵的阶数，返回时存放特征值
// 2. CMatrix& mtxEigenVector - 返回时存放特征向量矩阵，其中第i列为与
// 数组dblEigenValue中第j个特征值对应的特征向量
// 3. double eps - 计算精度，默认值为0.000001
//
// 返回值：BOOL型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::JacobiEigenv2(double dblEigenValue[], CMatrix& mtxEigenVector, double eps /*= 0.000001*/)
{
int i,j,p,q,u,w,t,s;
double ff,fm,cn,sn,omega,x,y,d;

if (! mtxEigenVector.Init(m_nNumColumns, m_nNumColumns))
return FALSE;

for (i=0; i<=m_nNumColumns-1; i++)
{
mtxEigenVector.m_pData=1.0;
for (j=0; j<=m_nNumColumns-1; j++)
if (i!=j)
mtxEigenVector.m_pData=0.0;
}

ff=0.0;
for (i=1; i<=m_nNumColumns-1; i++)
{
for (j=0; j<=i-1; j++)
{
d=m_pData;
ff=ff+d*d;
}
}

ff=sqrt(2.0*ff);

Loop_0:

ff=ff/(1.0*m_nNumColumns);

Loop_1:

for (i=1; i<=m_nNumColumns-1; i++)
{
for (j=0; j<=i-1; j++)
{
d=fabs(m_pData);
if (d>ff)
{
p=i;
q=j;
goto Loop_2;
}
}
}

if (ff<eps)
{
for (i=0; i<m_nNumColumns; ++i)
dblEigenValue = GetElement(i,i);
return TRUE;
}

goto Loop_0;

Loop_2:

u=p*m_nNumColumns+q;
w=p*m_nNumColumns+p;
t=q*m_nNumColumns+p;
s=q*m_nNumColumns+q;
x=-m_pData;
y=(m_pData-m_pData)/2.0;
omega=x/sqrt(x*x+y*y);
if (y<0.0)
omega=-omega;

sn=1.0+sqrt(1.0-omega*omega);
sn=omega/sqrt(2.0*sn);
cn=sqrt(1.0-sn*sn);
fm=m_pData;
m_pData=fm*cn*cn+m_pData*sn*sn+m_pData*omega;
m_pData=fm*sn*sn+m_pData*cn*cn-m_pData*omega;
m_pData=0.0; m_pData=0.0;

for (j=0; j<=m_nNumColumns-1; j++)
{
if ((j!=p)&&(j!=q))
{
u=p*m_nNumColumns+j;
w=q*m_nNumColumns+j;
fm=m_pData;
m_pData=fm*cn+m_pData*sn;
m_pData=-fm*sn+m_pData*cn;
}
}

for (i=0; i<=m_nNumColumns-1; i++)
{
if ((i!=p)&&(i!=q))
{
u=i*m_nNumColumns+p;
w=i*m_nNumColumns+q;
fm=m_pData;
m_pData=fm*cn+m_pData*sn;
m_pData=-fm*sn+m_pData*cn;
}
}

for (i=0; i<=m_nNumColumns-1; i++)
{
u=i*m_nNumColumns+p;
w=i*m_nNumColumns+q;
fm=mtxEigenVector.m_pData;
mtxEigenVector.m_pData=fm*cn+mtxEigenVector.m_pData*sn;
mtxEigenVector.m_pData=-fm*sn+mtxEigenVector.m_pData*cn;
}

goto Loop_1;
}

[ 本帖最后由风花雪月于 2006-11-14 19:23 编辑 ]

forwardhah 发表于 2005-11-13 15:04

"赫申伯格矩阵全部特征值的QR方法"

没有求对应的特征向量啊

有没有求特征向量的代码？

[ 本帖最后由风花雪月于 2006-11-14 19:25 编辑 ]

aliao 发表于 2006-4-20 20:30

下面的几个函数没有定义或说明过它的作用： fabs()、GetElement() 不管怎么样还是要谢谢风花雪月……

baiyang 发表于 2006-6-16 15:40

谢谢风花雪月了

norman 发表于 2006-6-22 00:20

回复：（aliao）下面的几个函数没有定义或说明过它的作...

fabs()是得到浮点数的绝对值
是c/c++ 带的库函数

[ 本帖最后由风花雪月于 2006-8-24 03:30 编辑 ]

markailee1 发表于 2006-6-22 23:20

GetElement()

GetElement() 是MFC的类COleSafeArray 的成员，
void GetElement( long* rgIndices, void* pvData );

[ 本帖最后由风花雪月于 2006-11-14 19:25 编辑 ]

carcanet 发表于 2006-6-28 19:59

谢谢风花雪月

maintell 发表于 2006-6-28 23:56

好啊,测试下代码

幻觉发表于 2006-7-20 17:15

今晚试试吧

xinyuxf 发表于 2006-8-23 16:51

自己编好麻烦啊，还是用matlab方便些，呵呵
好东西，赞

renyusen 发表于 2006-9-3 22:04

我在matlab中用的jacobi算法，还是那个方便，不过要学学用这些语言来编了，佩服佩服

hbyyj2008 发表于 2006-11-14 08:16

谢谢楼主

风花雪月 发表于 2006-11-14 19:25

原帖由 forwardhah 于 2005-11-13 15:04 发表
"赫申伯格矩阵全部特征值的QR方法"

没有求对应的特征向量啊

有没有求特征向量的代码？

特征向量还是比较容易求得，自己弄一下吧

页: [1] 2

声振论坛's Archiver

[分享]求矩阵特征值和特征向量的一个小程序

回复：（风花雪月）[分享]求矩阵特征值和特征向量的一...

回复：（风花雪月）[分享]求矩阵特征值和特征向量的一...

回复：（aliao）下面的几个函数没有定义或说明过它的作...

GetElement()