函数:
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(t>0)
{
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 (d>0)
{
xy=1
if (b<0)
{
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 (Loop1<1)
{
return false
}
Loop1--
j=t+2
while (j<m)
{
A[j*n+j-2]=0
j++
}
j=t+3
while (j<m)
{
A[j*n+j-3]=0
j++
}
k=t
while (k<m-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
q=A[(t+1)*n+t]*(A[t*n+t]+A[(t+1)*n+t+1]-x)
r=A[(t+1)*n+t]*A[(t+2)*n+t+1]
}
if (p!=0 || q!=0 || r!=0)
{
if (p<0)
{
xy=-1
}
else
{
xy=1
}
s=xy*sqrt(p*p+q*q+r*r)
if (k!=t)
{
A[k*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=kj<mj++)
{
b=A[k*n+j]
c=A[(k+1)*n+j]
p=x*b+e*c
q=e*b+y*c
r=f*b+g*c
if (k!=m-2)
{
b=A[(k+2)*n+j]
p+=f*b
q+=g*b
r+=z*b
A[(k+2)*n+j]=r
}
A[(k+1)*n+j]=q
A[k*n+j]=p
}
j=k+3
if (j>m-2)
{
j=m-1
}
for (i=ti<j+1i++)
{
b=A[i*n+k]
c=A[i*n+k+1]
p=x*b+e*c
q=e*b+y*c
r=f*b+g*c
if (k!=m-2)
{
b=A[i*n+k+2]
p+=f*b
q+=g*b
r+=z*b
A[i*n+k+2]=r
}
A[i*n+k+1]=q
A[i*n+k]=p
}
}
k++
}
}
}
delete []A
return true
}
#include<stdio.h>int main()
{
int a[3][5],i,j,k,temp
//通过读取给3×5的数组赋值
for(i=0i<3i++)
for(j=0j<5j++)
scanf("%d",&a[i][j])
//一行一行的判断
for(k=0k<3k++)
{
//标准冒泡法
for (j = 0j <9j++)
{
for (i = 0i <9 - ji++)
{
if (a[k][i] <a[k][i + 1])
{
temp = a[k][i]
a[k][i] = a[k][i + 1]
a[k][i + 1] = temp
}
}
}
for(i=0i<3i++)
{
for(j=0j<5j++)
printf("%d\t",a[i][j])
printf("\n")
}
}