tezhenzhi.cpp

用幂法和反幂法计算矩阵特征值急相关问题的一实例程序

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

int max(int a,int b,int c)
{
	if(a>=b && a>=c)
		return a;
	else if(b>a && b>=c)
		return b;
	else 
		return c;
}

int max(int a,int b)
{
	if(a>=b)
		return a;
	else 
		return b;
}
 
int min(int a,int b)
{
	if(a>b)
		return b;
	else 
		return a;
} 

void main()             //主程序
{
	int i,j,t;
	int m=501;             
	double **A=new double*[5];
	double *u=new double[m];
	double *y=new double[m];
	double *s=new double[m];
	double bert1=0.0,bert2=0.0,bert3=0.0;
	double ent;
	double k;
	double *p=new double[39];
	double *ber=new double[39];
	int h;

/////////////////////////////////////////////////////////////////////////////////A阵的初始化
	for(i=0;i<5;i++)
	{
		*(A+i)=new double[m];
	}
	for(i=0;i<m;i++)
	{
		*(*A+i)=-0.064;
		*(*(A+1)+i)=0.16;
		*(*(A+2)+i)=(1.64-0.024*(i+1))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1));
		*(*(A+3)+i)=0.16;
    	*(*(A+4)+i)=-0.064;
	}

//////////////////////////////////////////////////////////求按模最大的特征值
	for(i=0;i<m;i++)
	{
		*(u+i)=1;
		*(y+i)=0;
	}
	do
	{
		ent=0;
		for(i=0;i<m;i++)
		{
			ent=ent+(*(u+i))*(*(u+i));
		}
		ent=sqrt(ent);

		for(i=0;i<m;i++)
		{
			*(y+i)=*(u+i)/ent;
		}

		for(i=0;i<m;i++)
		{
			*(u+i)=0;
			for(j=0;j<m;j++)
			{
				if(abs(i-j)<3)
				{
					*(u+i)=*(u+i)+(*(*(A+i-j+2)+j))*(*(y+j));
				}
			}
		}

		ent=bert1;

		bert1=0;
		for(i=0;i<m;i++) 
		{
			bert1=bert1+(*(y+i))*(*(u+i));
		}
	}
	while(fabs((bert1-ent)/bert1)>=0.0000000000001);

////////////////////////////////////////////////////////求另一个端点特征值
	for(i=0;i<m;i++)
	{
		*(*(A+2)+i)=*(*(A+2)+i)-bert1;
	}
	for(i=0;i<m;i++)
	{
		*(u+i)=1;
		*(y+i)=0;
	}
   
do
	{
		ent=0;
		for(i=0;i<m;i++)
		{
			ent=ent+(*(u+i))*(*(u+i));
		}
		ent=sqrt(ent);

		for(i=0;i<m;i++)
		{
			*(y+i)=*(u+i)/ent;
		}
		for(i=0;i<m;i++)
		{
			*(u+i)=0;
			for(j=0;j<m;j++)
			{
				if(abs(i-j)<3)
				{
					*(u+i)=*(u+i)+(*(*(A+i-j+2)+j))*(*(y+j));
				}
			}
		}
	
		ent=bert2;

		bert2=0;
		for(i=0;i<m;i++) 
		{
			bert2=bert2+(*(y+i))*(*(u+i));
		}
	}
	while(fabs((bert2-ent)/bert2)>=0.0000000000001);
	bert2=bert2+bert1;
	for(i=0;i<m;i++)
	{
		*(*(A+2)+i)=*(*(A+2)+i)+bert1;
	}
//////////////////////////////////////////////////////决定大小
	if((bert2-bert1)>0)
	{
		k=bert1;
		ent=bert2;
	}
	else
	{
		k=bert2;
		ent=bert1;
	}
	bert1=k;
	bert2=ent;
///////////////////////////////////////////////////////求离某值最近的特征值
	for(h=0;h<=39;h++)
	{
		for(i=0;i<m;i++)
		{
			*(*A+i)=-0.064;
			*(*(A+1)+i)=0.16;
			*(*(A+2)+i)=(1.64-0.024*(i+1))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1));
			*(*(A+3)+i)=0.16;
			*(*(A+4)+i)=-0.064;
		}

		*(ber+h)=0;
		*(p+h)=bert1+(h+1)*(bert2-bert1)/40;
		for(i=0;i<m;i++)
		{
			*(*(A+2)+i)=*(*(A+2)+i)-(*(p+h));
		}
		for(i=0;i<m;i++)
		{
			*(u+i)=1;
			*(y+i)=0;
		}
		for(i=1;i<=m;i++)            //Doolittle分解
		{
			for(j=i;j<=min((i+2),m);j++)
			{
				k=*(*(A+i-j+2)+j-1);
				for(t=max(1,i-2,j-2);t<=(i-1);t++)
				{
					k=k-(*(*(A+i-t+2)+t-1))*(*(*(A+t-j+2)+j-1));
				}
				*(*(A+i-j+2)+j-1)=k;
			}
			for(j=i+1;j<=min(i+2,m);j++)
			{
				k=*(*(A+j-i+2)+i-1);
				for(t=max(1,j-2,i-2);t<=(i-1);t++)
				{
					k=k-(*(*(A+j-t+2)+t-1))*(*(*(A+t-i+2)+i-1));
				}
				k=k/(*(*(A+2)+i-1));
				*(*(A+j-i+2)+i-1)=k;
			}
		}
		
		for(i=0;i<m;i++)
		{
			*(u+i)=1;
			*(y+i)=0;
		}
		
		do
		{
			ent=0;
			for(i=0;i<m;i++)
			{
				ent=ent+(*(u+i))*(*(u+i));
			}
			ent=sqrt(ent);
			
			for(i=0;i<m;i++)
			{
				*(y+i)=*(u+i)/ent;
			}
			
			
			*s=*y;
			for(i=2;i<=m;i++)
			{
				k=*(y+i-1);
				for(t=max(1,i-2);t<=(i-1);t++)
				{
					k=k-(*(*(A+i-t+2)+t-1))*(*(s+t-1));
				}
				*(s+i-1)=k;
			}
			*(u+m-1)=*(s+m-1)/(*(*(A+2)+m-1));
			for(i=m-1;i>=1;i--)
			{
				k=*(s+i-1);
				for(t=i+1;t<=min(i+2,m);t++)
				{
					k=k-(*(*(A+i-t+2)+t-1))*(*(u+t-1));
				}
				k=k/(*(*(A+2)+i-1));
				*(u+i-1)=k;
			}
			
			ent=*(ber+h);
			
			*(ber+h)=0;
			for(i=0;i<m;i++) 
			{
				*(ber+h)=*(ber+h)+(*(y+i))*(*(u+i));
			}
		}
		while(fabs((*(ber+h)-ent)/(*(ber+h)))>=0.0000000000001);
		*(ber+h)=1.0/(*(ber+h));
		*(ber+h)=*(ber+h)+(*(p+h));

    	for(i=0;i<m;i++)
		{
			*(*(A+2)+i)=*(*(A+2)+i)+(*(p+h));
		}
	}
	
	///////////////////////////////////////////////////////求按模最小特征值
	for(i=0;i<m;i++)
	{
		*(*A+i)=-0.064;
		*(*(A+1)+i)=0.16;
		*(*(A+2)+i)=(1.64-0.024*(i+1))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1));
		*(*(A+3)+i)=0.16;
		*(*(A+4)+i)=-0.064;
	}
    for(i=0;i<m;i++)
	{
		*(u+i)=1;
		*(y+i)=0;
	}
	for(i=1;i<=m;i++)            //Doolittle分解
	{
		for(j=i;j<=min((i+2),m);j++)
		{
			k=*(*(A+i-j+2)+j-1);
			for(t=max(1,i-2,j-2);t<=(i-1);t++)
			{
				k=k-(*(*(A+i-t+2)+t-1))*(*(*(A+t-j+2)+j-1));
			}
			*(*(A+i-j+2)+j-1)=k;
		}
		for(j=i+1;j<=min(i+2,m);j++)
		{
			k=*(*(A+j-i+2)+i-1);
			for(t=max(1,j-2,i-2);t<=(i-1);t++)
			{
				k=k-(*(*(A+j-t+2)+t-1))*(*(*(A+t-i+2)+i-1));
			}
			k=k/(*(*(A+2)+i-1));
			*(*(A+j-i+2)+i-1)=k;
		}
	}

	for(i=0;i<m;i++)
	{
		*(u+i)=1;
		*(y+i)=0;
	}
	
    do
	{
		ent=0;
		for(i=0;i<m;i++)
		{
			ent=ent+(*(u+i))*(*(u+i));
		}
		ent=sqrt(ent);

		for(i=0;i<m;i++)
		{
			*(y+i)=*(u+i)/ent;
		}


		*s=*y;
		for(i=2;i<=m;i++)
		{
			k=*(y+i-1);
			for(t=max(1,i-2);t<=(i-1);t++)
			{
				k=k-(*(*(A+i-t+2)+t-1))*(*(s+t-1));
			}
			*(s+i-1)=k;
		}
		*(u+m-1)=*(s+m-1)/(*(*(A+2)+m-1));
		for(i=m-1;i>=1;i--)
		{
			k=*(s+i-1);
			for(t=i+1;t<=min(i+2,m);t++)
			{
				k=k-(*(*(A+i-t+2)+t-1))*(*(u+t-1));
			}
			k=k/(*(*(A+2)+i-1));
			*(u+i-1)=k;
		}
	
		ent=bert3;

		bert3=0;
		for(i=0;i<m;i++) 
		{
			bert3=bert3+(*(y+i))*(*(u+i));
		}
	}
	while(fabs((bert3-ent)/bert3)>=0.0000000000001);
	bert3=1.0/bert3;
	printf("最小的特征值：%.12e\n",bert1);
	printf("最大的特征值：%.12e\n",bert2);
	printf("按模最小的特征值：%.12e\n",bert3);
	for(h=0;h<=39;h++)
	{		
		printf("离%.12e最近的特征值：%.12e\n",*(p+h),*(ber+h));	
	}
///////////////////////////////////////////////////求谱范数   
	if((fabs(bert2)-fabs(bert1))>0)
	{
		k=bert2;
	}
	else
	{
		k=bert1;
	}
	k=fabs(k/bert3);
	printf("A的（谱范数）条件数是：%.12e\n",k);
//////////////////////////////////////////////////求A的行列式
	ent=1;
	for(i=0;i<501;i++)
	{
		ent=ent*(*(*(A+2)+i)); 
	}
	printf("A的行列式为：%.12e\n",ent);
}