#include<stdio.h>
#include<math.h>
#define n 10
#define L 2000
#define eg 1.0e-12
//全局变量
double A[n][n],Mk[n][n],Q[n][n];
double re1,re2,im1,im2,s,t;
double lm[n][2]; //存储特征值的矩阵,第一列为实部,第二列为虚部
int num;//特征值的个数
//函数声明部分
void output(double A[][n]);
void trangle(double A[][n]);
void QR(int m);
void count_root(double a,double b,double c);
void lie_gauss(double lm);
//主函数
void main()
{
int i,j,k,m,p;
double a,b,c;
double RQ[n][n];
re1=re2=im1=im2=0;
m=n-1;
num=0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
if(i==j)
A[i][j]=1.5*cos((i+1)+1.2*(j+1));
else
A[i][j]=sin(0.5*(i+1)+0.2*(j+1));
}
printf("The original matrix is:\n");
output(A);
trangle(A);
printf("经过拟上三角化的A阵:\n");
output(A);
for(i=0;i<n;i++)
{
for(j=0;j<2;j++)
{
lm[i][j]=0; //存储特征值的矩阵清零
}
}
for(k=1;k<=L;k++)
{
start:if(fabs(A[m][m-1])<=eg)
{
lm[num++][0]=A[m][m];
m--;
if(m==0)
{
lm[num++][0]=A[0][0];
break;
}
if(m>0)
goto start;
if(m<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];
count_root(a,b,c);
if(m==1)
{
lm[num][0]=re1;
lm[num++][1]=im1;
lm[num][0]=re2;
lm[num++][1]=im2;
break;
}
if(fabs(A[m-1][m-2])<=eg)
{
count_root(a,b,c);
lm[num][0]=re1;
lm[num++][1]=im1;
lm[num][0]=re2;
lm[num++][1]=im2;
m-=2;
if(m==0)
{
lm[num++][0]=A[0][0];
break;
}
if(m>0)
goto start;
if(m<0)
break;
}
if(k==L)
{
printf("没有得到全部的特征值,得到了%d个特征值\n",num);
break;
}
else
{
s=-b;
t=c;
QR(m);
}
}
printf("对矩阵进行QR分解后得到的矩阵Q和R分别如下:\n");
printf("Q阵:\n");
output(Q);
printf("\n");
printf("R阵:\n");
output(Mk);
printf("\n");
for(i=0;i<n;i++)//两个矩阵相乘
{
for(j=0;j<n;j++)
{
RQ[i][j]=0;
for(p=0;p<=m;p++)
{
RQ[i][j]+=Mk[i][p]*Q[p][j];
}
}
}
printf("RQ阵:\n");
output(RQ);
printf("\n");
printf("特征值如下所示:\n");
for(i=0;i<n;i++)
{
if(lm[i][1]==0)
printf("%.12e\n",lm[i][0]);
else if(lm[i][1]>0)
printf("%.12e+%.12ei\n",lm[i][0],lm[i][1]);
else
printf("%.12e%.12ei\n",lm[i][0],lm[i][1]);
}
//output(A);
for(i=0;i<n;i++)
{
if(lm[i][1]==0)
{
printf("第%d个非零特征值所对应的特征向量如下:",i+1);
lie_gauss(lm[i][0]);
printf("\n");
}
}
printf("\n");
}
//打印方阵函数
void output(double A[][n])
{
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
printf("%5.3f ",A[i][j]);
if(j==(n-1))
printf("\n");
}
}
//对A阵拟上三角函数
void trangle(double A[][n])
{
int r,i,j;
double Dr,Cr,Hr,temp=0;
double Ur[n]={0},Vr[n]={0},Pr[n]={0},Qr[n]={0},Tr,Wr[n]={0};
for(r=0;r<n-2;r++)
{
Tr=0;
for(i=0;i<n;i++)
{
Pr[i]=0;
Qr[i]=0;
}
for(i=r+2;i<n;i++)
{
if(A[i][r]!=0)
break;
}
if(i==n)
continue;
else
{
//计算Dr
for(i=r+1;i<n;i++)
temp+=A[i][r]*A[i][r];
Dr=sqrt(temp);
temp=0;
//计算Cr
if(A[r+1][r]<=0)
Cr=Dr;
else
Cr=-Dr;
//计算Hr
Hr=Cr*Cr-Cr*A[r+1][r];
//计算向量Ur
for(i=0;i<r+1;i++)
Ur[i]=0;
Ur[r+1]=A[r+1][r]-Cr;
for(i=r+2;i<n;i++)
Ur[i]=A[i][r];
//计算向量Vr
for(i=0;i<n;i++)
Vr[i]=Ur[i]/Hr;
//计算向量Pr
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
Pr[i]+=A[j][i]*Vr[j];
}
//计算向量Qr
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
Qr[i]+=A[i][j]*Vr[j];
}
//计算Tr
for(i=0;i<n;i++)
Tr+=Pr[i]*Vr[i];
//计算Wr
for(i=0;i<n;i++)
Wr[i]=Qr[i]-Tr*Ur[i];
//新的A阵
for(i=0;i<n;i++)
for(j=0;j<n;j++)
A[i][j]=A[i][j]-Wr[i]*Ur[j]-Ur[i]*Pr[j];
}
}
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
if(fabs(A[i][j])<1.0e-12)
A[i][j]=0;
}
}
}
//带双步位移的QR分解法
void QR(int m)
{
int k,i,j,r,p;
double ur[n]={0},vr[n]={0},pr[n]={0},qr[n]={0},wr[n]={0},lr[n]={0},er[n]={0},zancun[n][n];
double dr,cr,hr,tr=0,temp=0;
for(i=0;i<=m;i++)
{
for(j=0;j<=m;j++)
{
if(i==j)
Q[i][j]=1;
else
Q[i][j]=0;
}
}
for(i=0;i<=m;i++)
for(j=0;j<=m;j++)
Mk[i][j]=0;
for(i=0;i<=m;i++)//两个矩阵相乘
{
for(j=0;j<=m;j++)
{
Mk[i][j]=0;
for(p=0;p<=m;p++)
{
Mk[i][j]+=A[i][p]*A[p][j];
}
}
}
for(i=0;i<=m;i++)
{
for(j=0;j<=m;j++)
Mk[i][j]=Mk[i][j]-s*A[i][j];
}
for(i=0;i<=m;i++)
Mk[i][i]=Mk[i][i]+t;
//对Mk作QR分解并求A(k+1)
for(r=0;r<m;r++)
{
tr=0;
for(i=0;i<n;i++)
{
vr[i]=0;
pr[i]=0;
qr[i]=0;
er[i]=0;
}
for(i=r+1;i<=m;i++)
{
if(Mk[i][r]!=0)
break;
}
if(i>m)
continue;
else
{
//计算dr
for(i=r;i<=m;i++)
temp+=Mk[i][r]*Mk[i][r];
dr=sqrt(temp);
temp=0;
//计算cr
if(Mk[r][r]<=0)
cr=dr;
else
cr=-dr;
//计算hr
hr=cr*cr-cr*Mk[r][r];
//向量ur
for(i=0;i<r;i++)
ur[i]=0;
ur[r]=Mk[r][r]-cr;
for(i=r+1;i<=m;i++)
ur[i]=Mk[i][r];
//向量lr
for(i=0;i<=m;i++)
lr[i]=ur[i]/hr;
//向量er
for(i=0;i<=m;i++)
{
for(j=0;j<=m;j++)
er[i]+=Q[i][j]*lr[j];
}
//计算Qr
for(i=0;i<=m;i++)
{
for(j=0;j<=m;j++)
zancun[i][j]=er[i]*ur[j];
}
for(i=0;i<=m;i++)
{
for(j=0;j<=m;j++)
Q[i][j]-=zancun[i][j];
}
//向量vr
for(i=0;i<=m;i++)
{
for(j=0;j<=m;j++)
vr[i]+=Mk[j][i]*lr[j];
}
//新的Mk阵
for(i=0;i<=m;i++)
for(j=0;j<=m;j++)
Mk[i][j]=Mk[i][j]-ur[i]*vr[j];
//向量pr
for(i=0;i<=m;i++)
for(j=0;j<=m;j++)
pr[i]+=A[j][i]*lr[j];
//向量qr
for(i=0;i<=m;i++)
for(j=0;j<=m;j++)
qr[i]+=A[i][j]*lr[j];
//计算tr
for(i=0;i<=m;i++)
tr+=pr[i]*lr[i];
//向量wr
for(i=0;i<=m;i++)
wr[i]=qr[i]-tr*ur[i];
//新的A阵
for(i=0;i<=m;i++)
for(j=0;j<=m;j++)
A[i][j]=A[i][j]-wr[i]*ur[j]-ur[i]*pr[j];
}
}
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
if(fabs(A[i][j])<1.0e-12)
A[i][j]=0;
}
}
}
//计算一元二次方程
void count_root(double a,double b,double c)
{
double t;
t=b*b-4*a*c;
if(t>=0)
{
re1=(-b+sqrt(t))/(2*a);
re2=(-b-sqrt(t))/(2*a);
im1=0;
im2=0;
}
else
{
re1=-b/(2*a);
re2=re1;
im1=sqrt(-t)/(2*a);
im2=-im1;
}
}
//列主元素Gauss消去法求解特征向量
void lie_gauss(double lm)
{
int i,j,k,ik;
double A[n][n],temp,m,x[n]={0},sum;
for(i=0;i<n;i++) //产生A-λ*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))-lm;
}
}
//消元过程
for(k=0;k<n-1;k++)
{
//选行号ik
ik=k;
for(i=k+1;i<n;i++)
{
if(fabs(A[i][k])>fabs(A[ik][k]))
ik=i;
}
//交换行
if(ik!=k)
{
for(j=k;j<n;j++)
{
temp=A[k][j]; //temp暂存中间变量
A[k][j]=A[ik][j];
A[ik][j]=temp;
}
}
//计算
for(i=k+1;i<n;i++)
{
m=A[i][k]/A[k][k];
for(j=k+1;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++)
printf("%.12e ",x[i]);
}
- 1
- 2
- 3
- 4
- 5
- 6
前往页