chap9_jacobi.h

jacobi 和QR 算法 用来求解特征值问题 自己编的 不知道是否好用

const int dimA1=5;

class eigenJacobi
{
public:
	eigenJacobi(int num);	//构造函数，传递函数为矩阵A的维数
    void no1();				//输入第一题矩阵
	void jacobi();			//用Jacobi方法求解
	void output();          //输出结果
		
	double abs(double ab)   //求一个数的绝对值
	{
		return ab>0?ab:-ab;
	}
	
private :

	double a[dimA1][dimA1];    //存储矩阵A，并且用于输出特征值
	double b[dimA1][dimA1];    //在void jacobi()函数中暂时存储新的矩阵A'
	double v[dimA1][dimA1];    //存储变换矩阵V，并且用于输出特征向量
	double v1[dimA1][2];      //在void jacobi()函数中暂时存储变换矩阵
		
	int n;                  //n为矩阵的维数    
	double e;		        //容许误差
};


eigenJacobi::eigenJacobi(int num)//构造函数，传递函数为矩阵A的维数
{
	n=num;	                //n为矩阵的维数，根据题目设置为4            
	e=0.0000001;            //容许误差设置为e-7
	
	//设置初始的V矩阵
	for(int i=1;i<=n;i++)
	{
		for (int j=1;j<=n;j++)
		{
			if(i==j)
				v[i][j]=1;
			else 
				v[i][j]=0;
		}
	}
}

void eigenJacobi::no1()//输入第一题矩阵
{
	a[1][1]=10;	a[1][2]=7;	a[1][3]=8;	a[1][4]=7;
	a[2][1]=7;	a[2][2]=5;	a[2][3]=6;	a[2][4]=5;
	a[3][1]=8;	a[3][2]=6;	a[3][3]=10;	a[3][4]=9;
	a[4][1]=7;	a[4][2]=5;	a[4][3]=9;	a[4][4]=10;
}


void eigenJacobi::jacobi()//用Jacobi方法求解
{
	double max=0;         //记录模最大的元素值
    int i1,j1;			  //记录模最大的元素所在行列值
	double fai;           //旋转角度fai

	do
	{
	//先求按模最大的非对角线元素
		max=0;
		for (int i=1;i<=n;i++)
		{
			for (int j=1;j<i;j++)
			{
				if(max<abs(a[i][j]))
				{
					max=abs(a[i][j]);
					i1=i;
					j1=j;
				} 
			}
		}
		
	
	//然后计算A矩阵的值：

		//计算fai的值：
		if(abs(a[i1][i1]-a[j1][j1])>e)
			fai=0.5*atan(2*a[i1][j1]/(a[i1][i1]-a[j1][j1]));
		else 
			fai=3.1415926535897/4.0;

	
		//计算A'矩阵各项元素，使用b数组进行暂时的存储
		b[i1][i1]=a[i1][i1]*pow(cos(fai),2)+a[j1][j1]*pow(sin(fai),2)
			+2*a[i1][j1]*cos(fai)*sin(fai);
		b[j1][j1]=a[j1][j1]*pow(cos(fai),2)+a[i1][i1]*pow(sin(fai),2)
			-2*a[i1][j1]*cos(fai)*sin(fai);
		
		for(int l=1;l<=n;l++)
		{
			if(l!=i1&&l!=j1)
			{
				b[l][i1]=b[i1][l]=a[i1][l]*cos(fai)+a[j1][l]*sin(fai);
				b[j1][l]=b[l][j1]=-a[i1][l]*sin(fai)+a[j1][l]*cos(fai);

				for(int m=1;m<=n;m++)
					if(m!=i1&&m!=j1)
					{
						b[m][l]=a[m][l];
					}
			}
			
		}

		b[j1][i1]=b[i1][j1]=0;
		
		//存储新的A矩阵(即A'矩阵，b数组中存储的元素)
		for (i=1;i<=n;i++)
		{
			for (int j=1;j<=n;j++)
			{
				a[i][j]=b[i][j];
			}
		}


  
	//最后累次相乘计算V的值，以得到对应的特征向量
	//由于V矩阵的特殊性（仅有i1和j1行列交叉的四个数非1或者0）
    //因此做矩阵乘法的时候可以只考虑与i1和j1相关的行列变动即可：

		for (i=1;i<=n;i++)
		{
			v1[i][0]=v[i][i1]*cos(fai)+v[i][j1]*sin(fai);
			v1[i][1]=v[i][j1]*cos(fai)-v[i][i1]*sin(fai);
		}

		for (i=1;i<=n;i++)
		{
			v[i][i1]=v1[i][0];
			v[i][j1]=v1[i][1];
		}


	}while(max>e);
	
}



void eigenJacobi::output()
{	
	//输出结果:
	for (int ii=1;ii<=n;ii++)
	{		
		cout<<"第"<<ii<<"个特征值为：  "<<a[ii][ii]<<endl;
		cout<<"对应的特征向量为：[";
		for(int jj=1;jj<=n;jj++)
			cout<<setw(11)<<v[jj][ii];
		cout<<"   ]'"<<endl;
		cout<<endl;
	}	
}