meller-幂法.txt

之前做的meller-幂法

#include<iostream>
#include<cmath>
#include<iomanip>
using namespace std;
#define N 3
//矩阵与向量相乘
/*void fun(float a[N][N],float u[N])
{
	for(int i=0;i<N;i++)
		for(int j=0;j<N;j++)
			u[i]=u[i]+a[i][j]*u[j];
}*/
void main()
{
	//初始化
	float a[N][N]={{1.0,1.0,0.5},{1.0,1.0,0.25},{0.5,0.25,2.0}};
	float u0[N]={1.0,1.0,1.0};
	float v1[N];
	float u1[N];
	

	double ep1,eps1,ep,eps,max_v,max_u,t(1.0);
	int n=0;
	cout<<"输入控制精度（ep,eps(max_v)）"<<endl;
	cin>>ep;
	cin>>eps;

	ep1=2*ep;

	cout<<"k"<<"      U(k)的规范化向量"<<"              max(V(k))"<<endl;
	cout<<n<<"      ";
	for(int k=0;k<N;k++)
		{
		cout<<setiosflags(ios_base::left)<<setw(10)<<u0[k];
		}
	cout<<endl;
	                      
	while(n<500 && (ep1>ep || eps1>eps) )
	{
	
	//	fun(a,u1);
		for(int k=0;k<N;k++)
		{
			v1[k]=0;
		}
		for(int i=0;i<N;i++)
		{	for(int j=0;j<N;j++)
				v1[i]=v1[i]+a[i][j]*u0[j];
		}
		
		n++;
		
		//求max(vk)
		max_v=v1[0];
		for( i=0;i<N-1;i++)
		{
			if( fabs(max_v) < fabs(v1[i+1]) )
				max_v=v1[i+1];
		}
		eps1=fabs(max_v-t);//求eps1
		t=max_v;//记录上次max(vk)

		//求uk
		for(int j=0;j<N;j++)
			u1[j]=v1[j]/max_v;
		//++++++ ||u(k)-u(k-1)||

		max_u=u1[0]-u0[0];
		for(j=0;j<N-1;j++)
		{
			if( fabs(max_u) < fabs(u1[j+1]-u0[j+1]) )
				max_u=fabs(u1[j+1]-u0[j+1]);
		}
		ep1=max_u;
		cout<<setw(4)<<n<<" ( ";
		for(k=0;k<N;k++)
		{
			u0[k]=u1[k];
			cout<<setiosflags(ios_base::left)<<setw(10)<<setprecision(6)<<u0[k];
		}
		cout<<")   ";
		cout<<setprecision(6)<<max_v<<endl;
	}
}