Question

c++如何求矩阵特征值

程序

Answer 1

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=k;j<m;j++)
{
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=t;i<j+1;i++)
{
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;
}

Answer 2

就是求特征值的话可以利用递归和高斯迭代法解高次方程，有了特征值之后代入，用Gauss_Seidel迭代法解线性方程组，即特征向量

Answer 3

//QR方法求矩阵特征值
/******
数据
-1 2 1
2 4 -1
1 1 -6
******/
#include<iostream>
#include<fstream>
#include<iomanip>
using namespace std;
#include<math.h>
#define N 3 //矩阵的维数
#define NUM 15 //QR分解次数
#define SNUM 13 //用来控制输出格式

void Print(long double A[N][N],long double Q[N][N],long double R[N][N]); //矩阵输出
void QR(long double A[N][N],long double Q[N][N],long double R[N][N]); //QR分解
void Multiplicate(long double A[N][N],long double R[N][N],long double Q[N][N]); //迭代，获得下一次矩阵A=QR

int main()
{
int i,j;
long double A[N][N];
long double Q[N][N]={0};
long double R[N][N]={0};

cout<<"*************************************************************"<<endl;
cout<<"输入初始矩阵:\n";
cout<<"-------------------------------------------------------------"<<endl;
for(i=0;i<N;i++)
for(j=0;j<N;j++)
cin>>A[i][j];
cout<<"*************************************************************"<<endl;

cout<<"输出每一步迭代过程: \n";
cout<<"*************************************************************"<<endl;
for(i=1;i<=NUM;i++)
{
QR(A,Q,R);
cout<<"第"<<i<<"步QR分解:\n";
cout<<"-------------------------------------------------------------"<<endl;
Print(A,Q,R);
Multiplicate(A,R,Q);
}

cout<<"矩阵特征值为:\n";
cout<<"-------------------------------------------------------------"<<endl;
for(i=0;i<N;i++) //输出特征值
cout<<"r["<<(i+1)<<"]="<<R[i][i]<<endl;
cout<<"*************************************************************"<<endl;
return 0;
}

void QR(long double A[N][N],long double Q[N][N],long double R[N][N])
{
int i,j,k,m;
long double temp;
long double a[N],b[N];
for(j=0;j<N;j++)
{
for(i=0;i<N;i++)
{
a[i]=A[i][j];
b[i]=A[i][j];
}
for(k=0;k<j;k++)
{
R[k][j]=0;
for(m=0;m<N;m++)
R[k][j]+=a[m]*Q[m][k];
for(m=0;m<N;m++)
b[m]-=R[k][j]*Q[m][k];
}
temp=0;
for(i=0;i<N;i++)
temp+=b[i]*b[i];
R[j][j]=sqrt(temp);
for(i=0;i<N;i++)
Q[i][j]=b[i]/sqrt(temp);
}
}

void Multiplicate(long double A[N][N],long double R[N][N],long double Q[N][N])
{
int i,j,k;
long double temp;
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{
temp=0;
for(k=0;k<N;k++)
temp+=R[i][k]*Q[k][j];
A[i][j]=temp;
}
}

void Print(long double A[N][N],long double Q[N][N],long double R[N][N])
{
int i,j;
cout<<left;
cout<<"矩阵A:\n";
for(i=0;i<N;i++){
for(j=0;j<N;j++)
cout<<setw(SNUM)<<A[i][j];
cout<<endl;
}

cout<<"-------------------------------------------------------------"<<endl;
cout<<"矩阵Q:\n";

for(i=0;i<N;i++){
for(j=0;j<N;j++)
cout<<setw(SNUM)<<Q[i][j];
cout<<endl;
}

cout<<"-------------------------------------------------------------"<<endl;
cout<<"矩阵R:\n";

for(i=0;i<N;i++){
for(j=0;j<N;j++)
cout<<setw(SNUM)<<R[i][j];
cout<<endl;
}
cout<<"*************************************************************"<<endl;
}

本文来自CSDN博客，转载请标明出处：http://blog.csdn.net/CONPERINT/archive/2011/02/10/6177919.aspx