20036127002姜永洪(反幂法).cpp

反密法在数值计算中的地位就不用我来一一列举了

#include<iostream.h>
#include<math.h>
#define  eps  1e-6
void LU(double**A,double**L,double**U,int n)
{
	int i,j,k;
	for(k=0;k<n;k++)  //LU分解    
	{
		L[k][k]=1;
		for(j=0;j<n;j++)
		{
			if(j>k) L[k][j]=0;
			if(j<k) U[k][j]=0;
		}
		for(j=k;j<n;j++)
		{
			for(i=0;i<k;i++)
			    A[k][j]-=L[k][i]*U[i][j];			
			U[k][j]=A[k][j];
		}
		for(i=k+1;i<n;i++)
		{
			for(j=0;j<k;j++)
				A[i][k]-=L[i][j]*U[j][k];
			L[i][k]=A[i][k]/U[k][k];
		}
	}
}
void solve(double**L,double**U,double*B,double*X,int n)
{
	int i,j,k;
	double*M;
	M=new double[n];
	for(k=0;k<n;k++)        // solve  X
	{
		M[k]=B[k];
		for(j=0;j<k;j++)
		{
			M[k]-=L[k][j]*M[j];
		}
	}
	for(k=n-1;k>-1;k--)
	{
		for(i=k+1;i<n;i++)
			M[k]-=U[k][i]*X[i];
		X[k]=M[k]/U[k][k];
	}
	delete[]M;
}
void input(double**A,int n)
{
	int i,j;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			cout<<"A["<<i+1<<"]["<<j+1<<"]=";
		    cin>>A[i][j];
		}
		cout<<endl;
	}
}
void main()
{
	int i,n,k;
	double**A,**L,**U,*B,*X,*Y,*U1,*U2,t,d,T,T1,tt,temp;
	cout<<"input the index of matrixa"<<endl;
	cout<<"index=";
	cin>>n;
	A=new double*[n];
	U=new double*[n];
	L=new double*[n];
	U1=new double[n];
	U2=new double[n];	
	B=new double[n];
	X=new double[n];
	Y=new double[n];
  	for(i=0;i<n;i++)
	{
		A[i]=new double[n];
		L[i]=new double[n];
		U[i]=new double[n];
	}
//	input(A,n);
	A[0][0]=1;A[0][1]=-1;A[0][2]=2;
	A[1][0]=-2;A[1][1]=0;A[1][2]=5;
	A[2][0]=6;A[2][1]=-3;A[2][2]=6;
	LU(A,L,U,n);
	for(i=0;i<n;i++)
	{
		U1[i]=1;
	}
    T=0;
	do
	{   
		t=0;
		tt=0;
		for(i=0;i<n;i++)         //确定 | A_r| =max| A_i|
		{	
			if(fabs(U1[i])>t)
			{
				t=fabs(U1[i]);
				temp=U1[i];
			}
		}
		if(temp>0)               //t_k=sgn(A_r)/(U_k)_r
			k=1;
		if(temp<0)
			k=-1;
		for(i=0;i<n;i++)
		{
			Y[i]=U1[i]/t;        //Y_k-1=U_k-1/|A_r|
		}
		solve(L,U,Y,U2,n);
		for(i=0;i<n;i++)
		{
			if(fabs(U2[i])>fabs(tt))
			{
				tt=U2[i];
			}
			U1[i]=U2[i];
		}
		T=k/tt;
		d=fabs(T-T1);
		T1=T;
		cout<<"................"<<T<<endl;
	}while(d>eps);
	cout<<"Output :"<<endl;
	cout<<" 特征值："<<endl;
	cout<<"T="<<T<<endl;
	cout<<" 特征向量："<<endl;
	for(i=0;i<n;i++)
	{
//		cout<<"Y["<<i+1<<"]="<<Y[i]<<"     ";
        cout<<Y[i]<<"    ";
	}
	cout<<endl;
	delete[]A;delete[]B;delete[]L;
	delete[]Y;delete[]X;delete[]U;
}