#include<iostream>
#include<math.h>
#include<iomanip>
using namespace std;
#define N 10
#define L 2000//迭代次数
double A[N][N];
double re1,re2,im1,im2,d_s,d_t;
double tzz[N][2];//特征值矩阵，第1列为实部，第2列为虚部
int tzz_num;//特征值个数
void triangle(double B[][N]);//拟上三角化
void QR(int m);//QR分解M矩阵并求出相应的Ak
void root(double a,double b,double c);//求一元二次方程
void gauss(double value);//列主元高斯消去法求解特征向量
void triangle(double B[][N])//拟上三角化
{
int i,j,r;
double sum_temp,d,c,h,t,div_uh[N],u[N],p[N],q[N],w[N];
for(r=0;r<N-2;r++)
{
sum_temp=0;//初始化sum_temp,t及p,q
t=0;
for(i=0;i<N;i++)
{
p[i]=0;
q[i]=0;
}
for(i=r+2;i<N;i++)
{
sum_temp+=B[i][r]*B[i][r];
}
if(sum_temp==0)
continue;
else
{
sum_temp+=B[r+1][r]*B[r+1][r];
d=sqrt(sum_temp);//求d
if(B[r+1][r]>0)//求c
c=-d;
else
c=d;
h=c*c-c*B[r+1][r];//求h
for(i=0;i<r+1;i++)//求u
u[i]=0;
u[r+1]=B[r+1][r]-c;
for(i=r+2;i<N;i++)
u[i]=B[i][r];
for(i=0;i<N;i++)//求div_uh
div_uh[i]=u[i]/h;
for(i=0;i<N;i++)//求p
{
for(j=0;j<N;j++)
p[i]=p[i]+B[j][i]*div_uh[j];//注意B的转置
}
for(i=0;i<N;i++)//求q
{
for(j=0;j<N;j++)
q[i]=q[i]+B[i][j]*div_uh[j];
}
for(i=0;i<N;i++)//求t
t+=p[i]*div_uh[i];
for(i=0;i<N;i++)//求w
w[i]=q[i]-t*u[i];
for(i=0;i<N;i++)//求A(r+1)
{
for(j=0;j<N;j++)
{
B[i][j]=B[i][j]-w[i]*u[j]-u[i]*p[j];
}
}
}
}
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
if(fabs(B[i][j])<1.0e-12)
B[i][j]=0;
}
}
printf("经过拟上三角化的矩阵如下：\n");
for(i=0;i<N;i++)
{
printf("\n");
for(j=0;j<N;j++)
{
printf("%f",B[i][j]);
printf(" ");
}
}
cout<<endl;
}
void root(double a,double b,double c)//求解一元二次方程
{
if((b*b-4*a*c)>=0)
{
re1=(-b+sqrt(b*b-4*a*c))/(2*a);
re2=(-b-sqrt(b*b-4*a*c))/(2*a);
im1=0;
im2=0;
}
if((b*b-4*a*c)<0)
{
re1=-b/(2*a);
re2=re1;
im1=sqrt(4*a*c-b*b)/(2*a);
im2=-sqrt(4*a*c-b*b)/(2*a);
}
}
void QR(int m)//QR分解法
{
int i,j,r,k;
double c,h,d,sum_temp,temp,w[N],u[N],p[N],div_uh[N],M[N][N],B[N][N],temp_A[N][N],Q[N][N],tran_Q[N][N];
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
if(i==j)
Q[i][j]=1;
else
Q[i][j]=0;
}
}
for(i=0;i<N;i++)
for(j=0;j<N;j++)
temp_A[i][j]=0;
for(i=0;i<N;i++)
for(j=0;j<N;j++)
for(k=0;k<N;k++)
temp_A[i][j]+=A[i][k]*A[k][j];
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
if(i<=m&&j<=m)
{
if(i==j)
M[i][j]=temp_A[i][j]-d_s*A[i][j]+d_t;
else
M[i][j]=temp_A[i][j]-d_s*A[i][j];
}
else
M[i][j]=0;
}
}
for(r=0;r<N-1;r++)
{
for(k=r+1;k<N;k++)
{
if (A[k][r]!=0)
break;
}
if (k==N)
{
continue;
}
sum_temp=0;
for (i=r;i<N;i++)
{
sum_temp+=M[i][r]*M[i][r];
}
d=sqrt(sum_temp);
if (M[r][r]>0)
c=-d;
else
c=d;
h=c*c-c*M[r][r];
for(i=0;i<r;i++)
{
u[i]=0;
}
u[r]=M[r][r]-c;
for(i=r+1;i<N;i++)
u[i]=M[i][r];
for(i=0;i<N;i++)
p[i]=0;
for(i=0;i<N;i++)
div_uh[i]=u[i]/h;
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
p[i]+=M[j][i]*div_uh[i];
}
}
for(i=0;i<N;i++)
w[i]=0;
for(i=0;i<N;i++)
for(j=0;j<N;j++)
w[i]+=Q[i][j]*u[j];
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
Q[i][j]=Q[i][j]-w[i]*div_uh[j];
M[i][j]=M[i][j]-u[i]*p[j];
}
}
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
tran_Q[i][j]=Q[j][i];
}
}
}
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{
temp=0;
for(k=0;k<N;k++)
{
temp+=tran_Q[i][k]*A[k][j];
}
B[i][j]=temp;
}
for(i=0;i<N;i++)
for(j=0;j<N;j++)
A[i][j]=B[i][j];
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{
temp=0;
for(k=0;k<N;k++)
{
temp+=A[i][k]*Q[k][j];
}
B[i][j]=temp;
}
for(i=0;i<N;i++)
for(j=0;j<N;j++)
A[i][j]=B[i][j];
}
void gauss(double value)//高斯方法求解特征向量
{
int i,j,k,ik;
double m,a[N][N],x[N],sum,temp;
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
if(i!=j)
a[i][j]=sin(0.5*(i+1)+0.2*(j+1));
else
a[i][j]=1.5*cos((i+1)+1.2*(j+1))-value;
}
}
for(k=0;k<N-1;k++)
{
ik=k;
for(i=k+1;i<N;i++)
{
if(fabs(a[ik][k])<fabs(a[i][k]))
ik=i;
}
if(ik!=k)
for(j=k+1;j<N;j++)
{
temp=a[k][j];
a[k][j]=a[ik][j];
a[ik][j]=temp;
}
}
for(k=0;k<N-1;k++)
{
for(i=k+1;i<N;i++)
{
m=a[i][k]/a[k][k];
for(j=k;j<N;j++)
a[i][j]=a[i][j]-m*a[k][j];
}
}
x[N-1]=1;
for(k=N-2;k>=0;k--)
{
sum=0;
for(j=k+1;j<N;j++)
{
sum+=a[k][j]*x[j];
}
x[k]=-sum/a[k][k];
}
for(i=0;i<N;i++)
{
cout<<x[i]<<" ";
}
}
void main()
{
int i,j,k,m;
double a,b,c;
re1=0;
re2=0;
im1=0;
im2=0;
m=N-1;
tzz_num=0;
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
if(i!=j)
A[i][j]=sin(0.5*(i+1)+0.2*(j+1));
else
A[i][j]=1.5*cos((i+1)+1.2*(j+1));
}
}
cout<<"原始的A矩阵为"<<endl;//输出原始A矩阵
for(i=0;i<N;i++)
{
printf("\n");
for(j=0;j<N;j++)
{
printf("%f",A[i][j]);
printf(" ");
}
}
cout<<endl;
triangle(A);
for(i=0;i<N;i++)
{
for(j=0;j<2;j++)
{
tzz[i][j]=0;
}
}
for(k=1;k<=L;k++)
{
if(fabs(A[m][m-1])<=1.0e-12)
{
tzz[tzz_num++][0]=A[m][m];
m--;
}
if(m==0)
{
tzz[tzz_num++][0]=A[0][0];
break;
}
a=1;
b=-(A[m-1][m-1]+A[m][m]);
c=A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1];
if(m==1)
{
root(a,b,c);
tzz[tzz_num][0]=re1;
tzz[tzz_num++][1]=im1;
tzz[tzz_num][0]=re2;
tzz[tzz_num++][1]=im2;
break;
}
if(fabs(A[m-1][m-2])<=1.0e-12)
{
root(a,b,c);
tzz[tzz_num][0]=re1;
tzz[tzz_num++][1]=im1;
tzz[tzz_num][0]=re2;
tzz[tzz_num++][1]=im2;
m=m-2;
continue;
}
if(k>L)
{
cout<<"未得到全部特征值，得到了"<<tzz_num<<"个特征值"<<endl;
break;
}
else
{
d_s=-b;
d_t=c;
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
if((i>=m+1)||(j>=m+1))
{
A[i][j]=0;
}
}
}
}
QR(m);
continue;
}
printf("特征值如下\n");
for(i=0;i<N;i++)
{
printf("\n");
if(tzz[i][1]==0)
{
cout<<tzz[i][0];
}
if(tzz[i][1]>0)
{
cout<<tzz[i][0]<<"+"<<tzz[i][1]<<"i";
}
if(tzz[i][1]<0)
{
cout<<tzz[i][0]<<tzz[i][1]<<"i";
}
}
printf("\n");
printf("特征向量如下\n");
for(i=0;i<N;i++)
{
printf("\n");
if(tzz[i][1]==0)
{
gauss(tzz[i][0]);
}
}
printf("\n");
system("pause");
}
- 1
- 2
前往页