c语言求矩阵特征向量函数 c语言求矩阵的迹( 五 )


a[u]=fm*cn+a[w]*sn;
a[w]=-fm*sn+a[w]*cn;
}
}
for (i=0; i=n-1; i++)
{
u=i*n+p;
w=i*n+q;
fm=v[u];
v[u]=fm*cn+v[w]*sn;
v[w]=-fm*sn+v[w]*cn;
}
}
return(1);
}
Jacobi法求实对称矩阵的特征值与特征向量(c语言 , 不要复制的,能注释清楚更好)void Jacobi(int n,float (*a)[7][7],float (*s)[8][8]) //n为矩阵阶数,a为输入矩阵c语言求矩阵特征向量函数,s为输出矩阵{ int i,j,i1,l,iq,iql,ip; float g,s1,s2,s3,v1,v2,v3,u,st,ct;for(i=0;i=n;i++) {for(j=0;j=i;j++){if((i-j)==0)(*s)[i][j]=1;else{(*s)[i][j]=0.0;(*s)[j][i]=0.0;}} } g=0; for(i=1;in;i++) {i1=i-1;for(j=0;j=i1;j++)g=g+2*(*a)[i][j]*(*a)[i][j]; } s1=sqrt(g); s2=min/(n)*s1; s3=s1;do {s3=s3/(n);do{l=0;for(iq=1;iqn;iq++){iql=iq-1;for(ip=0;ip=iql;ip++){if(fabs((*a)[ip][iq])=s3){l=1;v1=(*a)[ip][ip];v2=(*a)[ip][iq];v3=(*a)[iq][iq];u=0.5*(v1-v3);if(u==0)g=1;if(fabs(u)=1e-10)g=-(u/fabs(u)*l)*v2/sqrt(v2*v2+u*u);st=g/sqrt(2*(l+sqrt(l-g*g)));ct=sqrt(1-st*st);for(i=0;in;i++){g=(*a)[i][ip]*ct-(*a)[i][iq]*st;(*a)[i][iq]=(*a)[i][ip]*st+(*a)[i][iq]*ct;(*a)[i][ip]=g;g=(*s)[i][ip]*ct-(*s)[i][iq]*st;(*s)[i][iq]=(*s)[i][ip]*st+(*s)[i][iq]*ct;(*s)[i][ip]=g;}for(i=0;in;i++){(*a)[ip][i]=(*a)[i][ip];(*a)[iq][i]=(*a)[i][iq];}g=2*v2*st*ct;(*a)[ip][ip]=v1*ct*ct+v3*st*st-g;(*a)[iq][iq]=v1*st*st+v3*ct*ct+g;(*a)[ip][iq]=(v1-v3)*st*ct+v2*(ct*ct-st*st);(*a)[iq][ip]=(*a)[ip][iq];}}}}while(l==1); }while(s3s2); }
c++如何求矩阵特征值c++求矩阵特征值方法:
函数:
Matrix_EigenValue(double *K1,int n,int LoopNumber,double Error1,double *Ret)
K1:n阶方阵
n:方阵K1的阶数
LoopNumber:在误差无法保证能得到结果时运算的最大次数
Error1:误差控制变量
Ret:返回的一个n*2的矩阵 。矩阵的每一行是求得的特征值,第一列代表特征值实数部分,第二列代表虚数部分
函数成功返回True,失败返回False
特别说明:
Matrix_Hessenberg:把n阶方阵K1化为上三角Hessenberg矩阵,其中A储存上三角Hessenberg矩阵
源代码:
bool Matrix_EigenValue(double *K1,int n,int LoopNumber,double Error1,double *Ret)
{
int i,j,k,t,m,Loop1;
double b,c,d,g,xy,p,q,r,x,s,e,f,z,y,temp,*A;
A=new double[n*n];
Matrix_Hessenberg(K1,n,A);
m=n;
Loop1=LoopNumber;
while(m!=0)
{
t=m-1;
while(t0)
{
temp=abs(A[(t-1)*n+t-1]);
temp+=abs(A[t*n+t]);
temp=temp*Error1;
if (abs(A[t*n+t-1])temp)
{
t--;
}
else
{
break;
}
}
if (t==m-1)
{
Ret[(m-1)*2]=A[(m-1)*n+m-1];
Ret[(m-1)*2+1]=0;
m-=1;
Loop1=LoopNumber;
}
else if(t==m-2)
{
b=-A[(m-1)*n+m-1]-A[(m-2)*n+m-2];
c=A[(m-1)*n+m-1]*A[(m-2)*n+m-2]-A[(m-1)*n+m-2]*A[(m-2)*n+m-1];
d=b*b-4*c;
y=sqrt(abs(d));
if (d0)
{
xy=1;
if (b0)
{
xy=-1;
}
Ret[(m-1)*2]=-(b+xy*y)/2;
Ret[(m-1)*2+1]=0;
Ret[(m-2)*2]=c/Ret[(m-1)*2];
Ret[(m-2)*2+1]=0;
}
else
{
Ret[(m-1)*2]=-b/2;
Ret[(m-2)*2]=-b/2;
Ret[(m-1)*2+1]=y/2;
Ret[(m-2)*2+1]=-y/2;
}
m-=2;
Loop1=LoopNumber;
}
else
{
if (Loop11)
{
return false;
}
Loop1--;
j=t+2;
while (jm)
{
A[j*n+j-2]=0;
j++;
}
j=t+3;
while (jm)
{
A[j*n+j-3]=0;
j++;
}
k=t;
while (km-1)
{
if (k!=t)
{
p=A[k*n+k-1];
q=A[(k+1)*n+k-1];
if (k!=m-2)
{
r=A[(k+2)*n+k-1];
}
else
{
r=0;
}
}
else
{
b=A[(m-1)*n+m-1];
c=A[(m-2)*n+m-2];
x=b+c;
y=b*c-A[(m-2)*n+m-1]*A[(m-1)*n+m-2];
p=A[t*n+t]*(A[t*n+t]-x)+A[t*n+t+1]*A[(t+1)*n+t]+y;

推荐阅读