mifa.cpp

这是数值分析中幂法的实现程序

/***********************************************************/
/*     幂法求矩阵的最大特征值和特征向量Mifa1(),Mifa2()     */
/*			反幂法求矩阵的最小特征值和特征向量Fmifa()     */
/*	
/*
/***********************************************************/
#include "stdio.h"
#include "math.h"
#define n 3
double x[n];//方程组的解
////////////打印函数/////////////
void print(double A[n][n])
{
	printf("A is:\n");
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
		{	
			printf("%8.3f,",A[i][j]);
			if(j==n-1) printf("\n");
		}
}
int sign(double x)
{
	if(x>0) return 1;
	else if(x<0) return -1;
	else return 0;
}
void print1(char ar[2],double *arr)
{
	printf("%s is:",ar);
	for(int j=0;j<n;j++)				
	{	
		if(j==0) printf("\t[");
		printf("%14.9f",*(arr+j));
		if(j!=n-1) printf(",");
		else printf("\t]");
	}
	printf("\n");
}
///////求数组的最大值////////////
int Max2(int k1,double *array)
{
	double max=0;
	int maxnum=0;
	for(int kk=k1-1;kk<n;kk++)
		if(max<fabs(*(array+kk))) {max=fabs(*(array+kk));maxnum=kk;}
	printf("maxnum is %d,max is %5.3f\n",maxnum,max);
	return maxnum;
}
//////////////doolittle分解算法//////////////
void Doolittle(double A[n][n],double b[n])
{
	int M[n];//用于记录第k个主元素所在的行号
	int k;
	int i,t,j;
	double S[n];//中间量
	double y[n];
	double temp=0;
	double B[n][n],b1[n];
	int ik;
	for(i=0;i<n;i++)
		{
			for(j=0;j<n;j++)
				B[i][j]=A[i][j];
			b1[i]=b[i];
		}
	for(k=1;k<=n;k++)
	{
		//计算中间量si
		for(i=k;i<=n;i++)
		{
			for(t=1,temp=0;t<=k-1;t++)
			{
				temp=temp+A[i-1][t-1]*A[t-1][k-1];
			}
			S[i-1]=A[i-1][k-1]-temp;
			A[i-1][k-1]=S[i-1];
			//print();
		}
		ik=Max2(k,S);//(2)选行号ik
		M[k-1]=ik;
		printf("ik=%d,k-1=%d\n",ik,k-1);
		if(ik!=k-1)//(3)
		{
			for(int j=1;j<=n;j++)
			{temp=A[k-1][j-1];A[k-1][j-1]=A[ik][j-1];A[ik][j-1]=temp;}
			//print();	
		}
		//////(4)计算ukk,ukj,lik//////
		
		for(int j=k+1;j<=n;j++)
		{
			for(t=1,temp=0;t<=k-1;t++)
			{
				temp=temp+A[k-1][t-1]*A[t-1][j-1];
			}
			A[k-1][j-1]=A[k-1][j-1]-temp;//计算ukj
			A[j-1][k-1]=A[j-1][k-1]/A[k-1][k-1];//计算lik
			//print();
		}
		
	}
	//print(A);
//	print2("M",M);
	//求Qb
	for(k=1;k<=n;k++)
	{
		int tt;
		tt=M[k-1];
		temp=b[k-1];b[k-1]=b[tt];b[tt]=temp;
	}
	print1("b",b);
	//求解Ly=Qb和Ux=y
	y[0]=b[0];//y1
	for(i=2;i<=n;i++)
	{
		for(t=1,temp=0;t<=i-1;t++)
		{
			temp=temp+A[i-1][t-1]*y[t-1];
		}
		y[i-1]=b[i-1]-temp;//yi
	}
	print1("y",y);
	x[n-1]=y[n-1]/A[n-1][n-1];//xn
	for(i=n-1;i>=1;i--)
	{
		for(t=i+1,temp=0;t<=n;t++)
		{
			temp=temp+A[i-1][t-1]*x[t-1];
		}
		x[i-1]=(y[i-1]-temp)/A[i-1][i-1];//xi
	}
	for(i=0;i<n;i++)
		{
			for(j=0;j<n;j++)
				A[i][j]=B[i][j];
			b[i]=b1[i];
		}
	print1("x",x);
}

///////求数组的最大值////////////
int arraymax(double *array)
{
	double max=0;
	int maxnum=0;
	for(int kk=0;kk<n;kk++)
		if(max<fabs(*(array+kk))) {max=fabs(*(array+kk));maxnum=kk;}
	printf("maxnum is %d,max is %8.6f\n",maxnum,max);
	return maxnum;
}
double neiji(double arr1[n],double arr2[n])
{
	double value=0;
	for(int i=0;i<n;i++)
		value=value+arr1[i]*arr2[i];
	return value;
}
//第一种迭代方法
void MiFa1(double A[n][n],double ee)
{
	int i,j;
	int count=0;
	double temp=0;
	double u[n];
	double y[n],yita,beita1,beita=0;
	for(i=0;i<n;i++)
	{
		if(i==0) u[i]=1;
		else u[i]=0;
	}
	print(A);
	//print();//print array G[n][n];
	while(1)
	{
			beita1=beita;
			yita=sqrt(neiji(u,u));
			for(i=0;i<n;i++)
				y[i]=u[i]/yita;//求yk-1
			for(i=0;i<n;i++)
			{
				temp=0;
				for(j=0;j<n;j++)	
					temp=temp+A[i][j]*y[j];
				u[i]=temp;//求uk
			}
			beita=neiji(y,u);
			print1("u",u);
			print1("y",y);
			printf("beita is %11.6f\n",beita);
			count++;
			if(fabs(beita-beita1)/fabs(beita)<=ee) break;
	}
	printf("count is:%d\n",count);
	printf("nanmuda is %11.6f\n",beita);
	print1("tezhenxiangliang",y);
}
//第二种迭代方法
void MiFa2(double A[n][n],double ee)
{
	int i,j;
	int count=0;
	double temp=0;
	int rnum;//模最大分量号
	double hr;//模最大分量
	double u[n];
	double y[n],beita1,beita=0;
	for(i=0;i<n;i++)
	{
		if(i==0) u[i]=1;
		else u[i]=0;
	}
	print(A);
	//print();//print array G[n][n];
	while(1)
	{
			beita1=beita;
			rnum=arraymax(u);
			hr=u[rnum];
			for(i=0;i<n;i++)
				y[i]=u[i]/fabs(hr);//求yk-1
			for(i=0;i<n;i++)
			{
				temp=0;
				for(j=0;j<n;j++)	
					temp=temp+A[i][j]*y[j];
				u[i]=temp;//求uk
			}
			beita=sign(hr)*u[rnum];
			print1("u",u);
			print1("y",y);
			printf("beita is %11.6f\n",beita);
			count++;
			if(fabs(beita-beita1)/fabs(beita)<=ee) break;						
	}
	printf("count is:%d\n",count);
	printf("nanmuda is %11.6f\n",beita);
	print1("tezhenxiangliang",y);
}

//反幂法
void Fmifa(double A[n][n],double ee)
{
	int i;
	int count=0;
	double temp=0;
	double u[n]={1,0,0};
	double y[n],yita,beita1,beita=0;
	//for(i=0;i<n;i++)
		//u[i]=1;
	print(A);
	//print();//print array G[n][n];
	while(1)
	{
			beita1=beita;
			yita=sqrt(neiji(u,u));
			//printf("yita is %12.6f\n",yita);
			for(i=0;i<n;i++)
				y[i]=u[i]/yita;//求yk-1
			//Gauss(A,y);//求uk
			Doolittle(A,y);
			print(A);
			//print1("x",x);
			for(i=0;i<n;i++)
				u[i]=x[i];
			beita=neiji(y,u);
			print1("u",u);
			print1("y",y);
			printf("beita is %11.6f\n",beita);
			count++;
			if(fabs(1/beita-1/beita1)/fabs(1/beita)<=ee) break;
	}
	printf("count is:%d\n",count);
	printf("nanmuda is %11.6f\n",1/beita);
	print1("tezhenxiangliang",y);

}
int main(int argc, char* argv[])
{
	double A[n][n]={6,-12,6,-21,-3,24,-12,-12,51};
	//double A1[n][n]={5,30,-48,3,14,-24,3,15,-25};
	double A2[n][n]={2,3,2,10,3,4,3,6,1};
	double A3[n][n]={-8,1,0,1,-8,-3,0,-3,-1};
	//double A[n][n]={{0.012,0.01,0.167},{1,0.8334,5.91},{3200,1200,4.2}};
	//double b[n]={0.6781,12.1,981};
	//MiFa1(A,0.0001);
	//MiFa1(A2,0.0001);
	Fmifa(A3,0.0001);
	//Fmifa(A2,0.0005);
	return 0;
}