如何用c语言写求矩阵的特征值和特征向量

Python027

如何用c语言写求矩阵的特征值和特征向量,第1张

方法1:推导出det(aA-I)=0的解析式,这应该是个四次方程,因为只有4阶,不是很困难的,写出后就可以用方程求根的方法求解(如newton迭代法)方法2:如果你是对角优势阵,也就是对角线上的值的绝对值,比同行所有其他元素的绝对值的和还大,可以通过局部旋转的方法把矩阵“能量”集中到对角线这个是方法,你可以自己去写一下试试~

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(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")

}

}