special.h

用幂法与反幂法求矩阵的最大特征值及最小特征值

#include<math.h>
#include<stdio.h>
#include <conio.h>
#include <process.h>
#define N 501

int Function(b,d,n,l,il,m)
int n,l,il,m;
double b[],d[];
{
	int ls,k,i,j,is,u,v;
	double p,t;
	if(il!=2*l+1)
		printf ("fail !\n");return (-2);
	ls=l;
	for (k=0;k<n-1;k++){
		p=0.0;
		for (i=k;i<=ls;i++){
			t=fabs (b[i*il]);
			if(t>p)  p=t;is=i;	
		}
		if(p+1.0==1.0)printf("fail!\n");return (0);
		for (j=0;j<m;j++) {
			u=k*m+j;v=is*m+j;
			t=d[u];
			d[u]=d[v];
			d[v]=t;
		}
		for (j=0;j<il;j++){
			u=k*il+j;
			v=is*il+j;
			t=b[u];
			b[u]=b[v];
			b[v]=t;
		}
		for (j=0;j<m;j++){
			u=k*m+j;
			d[u]/=b[k*il];
		}
		for (j=1;j<il;j++){
			u=k*il+j;
			b[u]/=b[k*il];
		}
	 	for (i=k+1;i<=ls;i++){
			t=b[i*il];
			for (j=0;j<m;j++){
				u=i*m+j;
				v=k*m+j;
				d[u]-=t*d[v];
			}
			for (j=1;j<il;j++){
				u=i*il+j;
				v=k*il+j;
				b[u-1]=b[u]-t*b[v];
			}
			u=i*il+il-1;
			b[u]=0.0;
			}
			if(ls!=n-1) ls++;
		}
		p=b[(n-1)*il];
		if (fabs(p)+1.0==1.0){
			printf("fail!\n");
			return(0);
		}
		for(j=0;j<m;j++){
			u=(n-1)*m+j;
			d[u]/=p;
		}
		ls=1;
		for (i=n-2;i>=0;i--){
			for (k=0;k<m;k++){
				u=i*m+k;
				for(j=1;j<=ls;j++){
					v=i*il+j;is=(i+j)*m+k;
					d[u]-=b[v]*d[is];
				}
			}
			if(ls!=il-1) ls++;
		}
	return (2);
}


void SolveProblem()
{
	FILE * fp;
	static double s[N][5],h[N][5],a[N], b[N],c[N],d,e,f,g,MAX,MIN;
	int n,m,i,j,k;
	if ((fp=fopen("File1.txt","w"))==NULL){
		printf("can't open file");
		exit(0);
	}
	//利用矩阵的特殊形态采用幂法求矩阵的最大特征植及其特征向量
	for (i=0;i<N;i++){
		a[i]=(20.2-0.3*i)*sin(0.2*i+0.2);
		b[i]=0.1;
	}
	for (j=0;j<2000;j++){
		d=f=0;
		for(k=0;k<N;k++)
			d+=b[k]*b[k];
		d=sqrt (d);
		for (m=0;m<N;m++)
			c[m]=b[m]/d;
		for (n=2;n<N-2;n++)
			b[n]=-c[n-2]+2*c[n-1]+c[n]*a[n]+2*c[n+1]-c[n+2];
		b[0]=a[0]*c[0]+2*c[1]-c[2];
		b[1]=2*c[0]+a[1]*c[1]+2*c[2]-c[3];
		b[N-2]=-c[N-4]+2*c[N-3]+a[N-2]*c[N-2]+2*c[N-1];
		b[N-1]=-c[N-3]+2*c[N-2]+a[N-1]*c[N-1];
		for (i=0;i<N;i++)
			f+=b[i]*c[i];
		if(fabs((f-g)/f)<1.0e-12){
			MAX=f;
			break;
		}
		g=f;
		}
	for(i=0;i<N;i++){
		if (c[i]<1.0e-12)
			c[i]=0;
	}
	fprintf (fp,"最大特征值 MAX=%e \n",f);
	fprintf (fp,"最大特征向量如下: \n",f);
	for (i=0;i<20;i++)  fprintf(fp,"Ymax[%2d]=%12e   Ymax[%2d]=%12e\n",i,c[i],i+480,c[480+i]);
//采用矩阵的特殊形态利用反幂法及矩阵位移求矩阵的最小特征值及其特征向量
	for(i=2;i<N;i++){
		s[i][2]=(20.2-0.3*i)*sin(0.2*i+0.2);
		s[i][1]=s[i][3]=2;
		s[i][0]=s[i][4]=-1;
	}
	s[0][3]=s[0][4]=s[1][4]=s[N-2][4]=s[N-1][3]=s[N-1][4]=0;
	s[0][0]=20.2*sin(0.2);
	s[1][1]=19.9*sin(0.4);
	s[1][0]=s[0][1]=s[1][2]=2;
	s[0][2]=s[1][3]=-1;
	for (i=0;i<N;i++)	b[i]=1;
	for (j=0;j<2000;j++) {
		d=f=0;
		for(k=0;k<N;k++)
			d+=b[k]*b[k];
		d=sqrt (d);
		for(m=0;m<N;m++) {
			c[m]=b[m]/d;
			b[m]=c[m];
			for(k=0;k<5;k++)
			h[m][k]=s[m][k];
		}
		Function(h,b,N,2,5,1);
		for (i=0;i<N;i++)  f+=b[i]*c[i];
		if(fabs((f-g)/f)<1.0e-12){
			MIN=1/f; 
			for(i=0;i<N;i++)
			if (fabs(b[i])<1.0e-12) b[i]=0;
			fprintf (fp,"最小特征值MIN=%e \n",1/f);
			fprintf (fp,"最小特征向量如下: \n",f);
			for (i=0;i<20;i++)
			fprintf(fp,"Ymin[%2d]=%12e   Ymin[%2d]=%12e\n",i,b[i],i+480,b[480+i]);
			break;
		}
		g=f ;
	}
//采用矩阵的特殊形态利用反幂法及矩阵的位移求矩阵的最接近MIN+(MAX-MIN)*i/40   i=1,2,3......39的特征值
	fprintf (fp,"前39个特征值如下: \n",f);
	for (n=1;n<40;n++){
		for(i=2;i<N;i++)  s[i][2]=(20.2-0.3*i)*sin(0.2*i+0.2)-MIN-(MAX-MIN)*n/40;
		s[0][0]=20.2*sin(0.2)-MIN-(MAX-MIN)*n/40;
		s[1][1]=19.9*sin(0.4)-MIN-(MAX-MIN)*n/40;
		for (i=0;i<N;i++)
		b[i]=1;
		for (j=0;j<2000;j++){
			d=f=0;
			for(k=0;k<N;k++)
				d+=b[k]*b[k];
			d=sqrt (d);
			for (m=0;m<N;m++){
				c[m]=b[m]/d;
				b[m]=c[m];
				for(k=0;k<5;k++)
				h[m][k]=s[m][k];
			}
			fprintf (fp,"前39个特征值如下: \n",f);
			Function(h,b,N,2,5,1);
			for (i=0;i<N;i++)   f+=b[i]*c[i];
			if(fabs((f-g)/f)<1.0e-12){
				fprintf (fp,"R[%2d]=%e \n",n,1/f+MIN+(MAX-MIN)*n/40);
				break;
			}
			g=f;
		}
	}
	fprintf(fp,"cond(A)2=%e",fabs(MAX/MIN));
	fclose(fp);
}