hemework.cpp

求解矩阵的最大最小特征值及对应的特征向量

#include <stdio.h>
#include <math.h>

int min(int x,int y,int z);
int max(int x,int y,int z);
double mod(double u[501],int n );							  //求向量的模
void matrixmultiply(double A[5][501],double y[501],double u[501]);//矩阵与向量相乘
double vectormultiply(double y[501],double u[501]);			  //求向量内积
void LUdecomposition(double c[5][501]);						  // LU分解
void Equationsolution(double c[5][501],double b[501],double x[501]);//求解方程组的解
double Maxeig(double c[5][501],double u[501],double y[501]);  //求按模最大的特征值
double Mineig(double c[5][501],double u[501],double y[501]);  //求按模最小的特征值

void main()
{
	int i,j;
	int t;
	double c[5][501];
	double A[5][501];
	double u[501];
	double y[501];
	double temp;
	double maximum;                                           //存放最大特征值
	double maximumsector[501];								  //最大特征值对应的特征向量
	double minimum;											  //存放最小特征值
	double minimumsector[501];								  //最小特征值对应的特征向量
	double modulemaximum;									  //存放按模最大的特征值
	double moduleminimum;									  //存放按模最小的特征值
	double approximation[39];
	double eigseries[39];									  //要求的39个特征值
	double cond;											  //条件数
	FILE* fp;
	fp=fopen("result.txt","w");								  //结果存放文件

	//////////////////////////////////////////矩阵及迭代向量初始化,将所给矩阵A的值存入c[5][501]中
	for(i=0;i<=4;i++)
	{
		for(j=0;j<=500;j++)
		{
			c[i][j]=0.;
		}
	}
	for(i=0;i<=500;i++)
	{
		c[2][i]=(20.5-0.3*(i+1))*sin(0.2*(i+1));
		u[i]=0.;
		y[i]=0.;
	}
	for(i=0;i<=498;i++)
	{
		c[0][i+2]=-1.;
		c[4][i]=-1.;
	}
	for(i=0;i<=499;i++)
	{
		c[1][i+1]=2.;
		c[3][i]=2.;
	}
	u[0]=y[0]=1.;

	//////////////////////////////////////////////////////////求最大最小特征值
	temp=Maxeig(c,u,y);
	for(i=0;i<=4;i++)
		{
			for(j=0;j<=500;j++)
			{
				if(i==2)
					A[i][j]=c[i][j]-temp;
				else
					A[i][j]=c[i][j];
			}
		}
	if(temp>0.)
	{
		maximum=temp;
		for(j=0;j<=500;j++)
		{
			maximumsector[j]=y[j];
		}
		
		for(j=0;j<=500;j++)
		{
			u[j]=0.;
			y[j]=0.;
		}
		u[0]=y[0]=1.;
		minimum=Maxeig(A,u,y)+temp;
		for(j=0;j<=500;j++)
		{
			minimumsector[j]=y[j];
		}
	}
	else
	{
		minimum=temp;
		for(j=0;j<=500;j++)
		{
			minimumsector[j]=y[j];
		}
		for(j=0;j<=500;j++)
		{
			u[j]=0.;
			y[j]=0.;
		}
		u[0]=y[0]=1.;
		maximum=Maxeig(A,u,y)+temp;
		for(j=0;j<=500;j++)
		{
			maximumsector[j]=y[j];
		}
	}
	modulemaximum=(fabs(maximum)>fabs(minimum))? maximum:minimum;
	fprintf(fp,"   最小特征值为%.12e\n",minimum);
	fprintf(fp,"   最小特征值相应的特征向量的前20个分量为:\n");
	for(j=0;j<=19;j++)
	{
		fprintf(fp,"   %.12e",minimumsector[j]);
		if((j+1)%4==0)
			fprintf(fp,"\n");
	}
	fprintf(fp,"   最小特征值相应的特征向量的后20个分量为:\n");
	for(j=19;j>=0;j--)
	{
		fprintf(fp,"   %.12e",minimumsector[500-j]);
		if((20-j)%4==0)
			fprintf(fp,"\n");
	}
	fprintf(fp,"   最大特征值为%.12e\n",maximum);
	fprintf(fp,"   最大特征值相应的特征向量的前20个分量为:\n");
	for(j=0;j<=19;j++)
	{
		fprintf(fp,"   %.12e",maximumsector[j]);
		if((j+1)%4==0)
			fprintf(fp,"\n");
	}
	fprintf(fp,"   最大特征值相应的特征向量的后20个分量为:\n");
	for(j=19;j>=0;j--)
	{
		fprintf(fp,"   %.12e",maximumsector[500-j]);
		if((20-j)%4==0)
			fprintf(fp,"\n");
	}
//	printf("maximum=%.12e	minimum=%.12e\n",maximum,minimum);
	//////////////////////////////////////////////////////////求按模最小的特征值以求条件数
	for(j=0;j<=500;j++)
		{
			u[j]=0.;
			y[j]=0.;
		}
	u[0]=y[0]=1.;
	moduleminimum=Mineig(c,u,y);
	for(i=0;i<=4;i++)
	{
		for(j=0;j<=500;j++)
		{
			c[i][j]=0.;
		}
	}
	for(i=0;i<=500;i++)
	{
		c[2][i]=(20.5-0.3*(i+1))*sin(0.2*(i+1));
	}
	for(i=0;i<=498;i++)
	{
		c[0][i+2]=-1.;
		c[4][i]=-1.;
	}
	for(i=0;i<=499;i++)
	{
		c[1][i+1]=2.;
		c[3][i]=2.;
	}
	printf("modulemimum=%.12e\n",moduleminimum);
	//////////////////////////////////////////////////////////求要求的39个特征值
	for(t=0;t<=38;t++)
	{
		approximation[t]=minimum+(t+1)/40.*(maximum-minimum);
	}
	for(t=0;t<=38;t++)
	{
		for(j=0;j<=500;j++)
		{
			u[j]=0.;
			y[j]=0.;
		}
		u[0]=y[0]=1.;
		for(i=0;i<=4;i++)
		{
			for(j=0;j<=500;j++)
			{
				if(i==2)
					A[i][j]=c[i][j]-approximation[t];
				else
					A[i][j]=c[i][j];
			}
		}
		eigseries[t]=Mineig(A,u,y)+approximation[t];
	}
	fprintf(fp,"   所要求的39个特征值为:\n");
	for(t=0;t<=38;t++)
	{
		fprintf(fp,"   第%d个要求的特征值为:%.12e\n",t+1,eigseries[t]);
	}
//		printf("%.12e\n",eigseries[t]);
	//////////////////////////////////////////////////////////求条件数
	cond=fabs(modulemaximum/moduleminimum);
	fprintf(fp,"   谱范数条件数为:  %.12e\n",cond);
	fclose(fp);
//	printf("cond=%.12e\n",cond);

}

void LUdecomposition(double c[5][501])
{
	int i,j,k,t;
	for(k=0;k<=500;k++)
	{
		for(j=k;j<=min(k+2,500,500);j++)
		{
			for(t=max(0,k-2,j-2);t<=k-1;t++)
			{
				c[k-j+2][j]-=c[k-t+2][t]*c[t-j+2][j];
			}
		}
		for(i=k+1;i<=min(k+2,500,500);i++)
		{
			for(t=max(0,i-2,k-2);t<=k-1;t++)
			{
				c[i-k+2][k]-=c[i-t+2][t]*c[t-k+2][k];
			}
			c[i-k+2][k]/=c[2][k];
		}
	}
}
void Equationsolution(double c[5][501],double b[501],double x[501])
{
	int i,t;
	for(i=1;i<=500;i++)
	{
		for(t=max(0,i-2,i-2);t<=i-1;t++)
		{
			b[i]-=c[i-t+2][t]*b[t];
		}
	}
	x[500]=b[500]/c[2][500];
	for(i=499;i>=0;i--)
	{
		x[i]=b[i];
		for(t=i+1;t<=min(i+2,500,500);t++)
		{
			x[i]-=c[i-t+2][t]*x[t];
		}
		x[i]/=c[2][i];	
	}
}
int min(int x,int y,int z)
{
	int min=x;
	if(min>=y)
		min=y;
	if(min>=z)
		min=z;
	return min;
}
int max(int x,int y,int z)
{
	int max=x;
	if(max<=y)
		max=y;
	if(max<=z)
		max=z;
	return max;
}
double mod(double u[501], int n)
{
	double sum=0;
	double result;
	for(int i=0;i<=(n-1);i++)
	{
		sum+=u[i]*u[i];
	}
	result=sqrt(sum);
	return result;
}
void matrixmultiply(double A[5][501],double y[501],double u[501])
{
	int i,t;
	for(i=0;i<=500;i++)
	{
		u[i]=0.;
	}
	for(i=0;i<=500;i++)
	{
		for(t=max(0,i-2,i-2);t<=min(500,i+2,i+2);t++)
		{
			u[i]+=A[i-t+2][t]*y[t];
		}
	}

}
double vectormultiply(double y[501],double u[501])
{
	double result=0.;
	int i;
	for(i=0;i<=500;i++)
	{
		result+=y[i]*u[i];
	}
	return result;
}
double Mineig(double c[5][501],double u[501],double y[501])
{
	int k,t;
	double result=1e20;
	double ita;
	double b[501];
	double beta[5001];
	LUdecomposition(c);
	for(k=1;k<=5000;k++)
	{
		ita=mod(u,501);
		for(t=0;t<=500;t++)
		{
			y[t]=u[t]/ita;
			b[t]=y[t];
		}
		Equationsolution(c,b,u);
		beta[k]=vectormultiply(y,u);
		if(k>=2&&(fabs(1./beta[k]-1./beta[k-1])/fabs(1./beta[k]))<1.e-12)
			break;
	}
	if(k<=5000)
	{
		result=1./beta[k];
	}
		return result;
}
double Maxeig(double c[5][501],double u[501],double y[501])
{
	int k,t;
	double result=1e20;
	double ita;
//	double b[501];
	double beta[5001];
	for(k=1;k<=5000;k++)
	{
		ita=mod(u,501);
		for(t=0;t<=500;t++)
		{
			y[t]=u[t]/ita;
//			b[t]=y[t];
		}
		matrixmultiply(c,y,u);
		beta[k]=vectormultiply(y,u);
		if(k>=2&&(fabs(beta[k]-beta[k-1])/fabs(beta[k]))<1.e-12)
			break;
	}
	if(k<=5000)
	{
		result=beta[k];
	}
		return result;
}