chap9_qr.h

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

const int dimA2=5;
const int max_step=100000;    //QR分解最大迭代次数

class eigenQR 
{
public:
	eigenQR(int num);		//构造函数
    void no2();				//输入第二题矩阵

	void qr();				//QR解法，其中调用了下面的void schmidt()，void exchange()函数
	
	void schmidt();			//施密特正交化，其中调用了下面的double multiple(int i,int k)
	                        //以及double norm2(int i)函数
    void exchange();        //交换相乘
    
	double compare();       //计算本次循环得到的矩阵与上次循环矩阵的差别
                            //利用元素之差的按模取最大值

	double multiple(int i,int k);//计算点乘
	double norm2(int i);    //计算二范数

	double abs(double ab)   //计算绝对值
	{
		return ab>0?ab:-ab;
	}
	
private :

	double a[dimA2][dimA2];    //存储矩阵A，并且用于输出特征值
	double b[dimA2][dimA2];    //在void schmidt()函数中暂时存储未单位化的向量b'
	
	double q[dimA2][dimA2];    //存储QR分解中的Q矩阵          
	double r[dimA2][dimA2];    //存储QR分解中的R矩阵
	
	double aa[dimA2][dimA2];   //用来存储QR分解之前的A矩阵元素
	double ee[dimA2][dimA2];   //用来存储每次Q矩阵的乘积（即书中的E矩阵）
	                           //并且用于最后输出矩阵的特征向量

	int n;                     //矩阵A的维数   
	double e;                  //容许误差
};


eigenQR::eigenQR(int num)//构造函数
{
	n=num;	                //n为矩阵的维数，根据题目设置为4            
	e=0.0000001;            //容许误差设置为e-7

	//设置初始的EE矩阵
	for (int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)		
		{
			if(i==j)
				ee[i][i]=1;
			else ee[i][j]=0;
		}	
	}	
}

void eigenQR::no2()//输入第二题矩阵
{
	a[1][1]=2;	a[1][2]=-1;	a[1][3]=0;	a[1][4]=0;
	a[2][1]=-1;	a[2][2]=2;	a[2][3]=-1;	a[2][4]=0;
	a[3][1]=0;	a[3][2]=-1;	a[3][3]=2;	a[3][4]=-1;
	a[4][1]=0;	a[4][2]=0;	a[4][3]=-1;	a[4][4]=2;

}


void eigenQR::qr()//QR解法
{

	//循环迭代max_step=100000次

	for (int ii=1;ii<=max_step;ii++)
	{
		schmidt();       //先进行施密特正交化，分解为Q和R矩阵
		exchange();      //然后交换相乘

		//如果两次循环求出的矩阵相差小于误差限的时候，停止循环，输出结果
		if(compare()<e)
			break;
	}

	//输出结果:
	cout<<"经过"<<ii<<"次迭代后，得到的结果如下所示:"<<endl<<endl;
	for (ii=1;ii<=n;ii++)
	{		
		cout<<"第"<<ii<<"个特征值为：  "<<a[ii][ii]<<endl;
		cout<<"对应的特征向量为：[";
		for(int jj=1;jj<=n;jj++)
			cout<<setw(11)<<ee[jj][ii];
		cout<<"   ]'"<<endl;
		cout<<endl;
	}	

}


void eigenQR::schmidt()//施密特正交化，分解为Q和R矩阵
{

//求Q矩阵：
	for (int i=1;i<=n;i++)
	{
		//计算b'向量中的各分量,用b数组进行存储
		for (int j=1;j<=n;j++)
		{
			b[j][i]=a[j][i];
			for (int k=1;k<=i-1;k++)
			{
				b[j][i]=b[j][i]-multiple(i,k)*q[j][k];
				//其中multiple(i,k)是a中第i列向量与q中第k列向量进行点乘，函数定义见后
			}
		}

		//计算Q中的元素，将b'向量单位化
		for (j=1;j<=n;j++)
		{
			q[j][i]=b[j][i]/norm2(i);
			//其中norm2(i)是求b'中第i列向量的二范数，函数定义见后
			
		}

	}


//求R矩阵：
	for (i=1;i<=n;i++)
	{
		//R对角线元素即为b'中第i列向量的二范数
		r[i][i]=norm2(i);
		
		//对角线右上的元素为a中第j列向量与q中第i列向量进行点乘
		for (int j=i+1;j<=n;j++)
		{
			r[i][j]=multiple(j,i);
		}

		//其余元素为0
		for (j=1;j<i;j++)
		{
			r[i][j]=0;
		}
	}
}

void eigenQR::exchange()//对QR进行交换相乘，并且累次相乘Q矩阵得到新的E矩阵
{

	//aa数组用来存储QR分解之前的A矩阵元素，用于之后compare函数中的比较
	//tt数组用来暂时存储上一次的E矩阵的值
	double tt[20][20];
	for (int i=1;i<=n;i++)
	{
		for (int j=1;j<=n;j++)
		{    
			aa[i][j]=a[i][j];
			tt[i][j]=ee[i][j];
		}
	}
	
	//计算新的A矩阵以及E矩阵
	for (i=1;i<=n;i++)
	{
		for (int j=1;j<=n;j++)
		{    
			a[i][j]=0;
			ee[i][j]=0;
			for (int k=1;k<=n;k++)
			{
				a[i][j]+=r[i][k]*q[k][j];      //RQ相乘
				ee[i][j]+=tt[i][k]*q[k][j];    //累次乘上新的Q矩阵
			}
		
		}
	}
}


double eigenQR::compare()//计算本次循环得到的矩阵与上次循环矩阵的差别
                         //利用元素之差的按模取最大值
{
	double temp=0;
	for(int i=1;i<=n;i++)
		for (int j=1;j<=n;j++)
			if ( temp<abs(a[i][j]-aa[i][j]) )
				temp=abs(a[i][j]-aa[i][j]);
	return temp;
}


double eigenQR::multiple(int i,int k)//a中第i列向量与q中第k列向量进行点乘
{   
	double temp=0;
	for (int j=1;j<=n;j++)
	{
		temp+=a[j][i]*q[j][k];
	}
	return temp; 
}



double eigenQR::norm2(int i)//计算b'中第i列向量的二范数
{   
	double temp=0;
	for (int j=1;j<=n;j++)
	{
		temp+=pow(b[j][i],2);
	}
	return sqrt(temp);
	
}