jacobi.cpp

雅克比迭代法求实对称矩阵的特征值和特征向量

#include "stdio.h"
#include "math.h"

bool main()
{
	int iNum = 3;
	int iMaxFN = 3;
	int nMaxIt = 200; //迭代次数
	double eps = 0.0001; //精度控制

	double EigenValue[4];
	double EigenVector[4][4];
	int i,j,p,q,u,w,t,s,l;
    double fm,cn,sn,omega,x,y,d;
	
	double iR1[3][3] = {
		{2,-1,3},
		{-1,2,-1},
		{3,-1,2}
	};
	
//	double iR1[4][4] = {
//		{2,-1,0,3},
//		{-1,2,-1,3},
//		{0,-1,2,3},
//		{3,3,3,3}
//	};
	l = 1;
	printf("\n");
    for (i=0; i <= iNum-1; i++)
    {
		EigenVector[i][i] = 1.0;
        for (j = 0; j <= iNum-1; j++)
			if (i != j) 
				EigenVector[i][j]=0.0;
    }

	printf("输入的矩阵：\n");
	for(i = 0;i < iNum;i++)
	{
		for(j = 0;j < iNum;j++)
			printf("%f ",iR1[i][j]);
		printf("\n");
	}
    
	while (l < nMaxIt)
	{ 
		fm=0.0;
        for (i = 1; i <= iNum - 1; i++)
		{
			for (j = 0; j <= i-1; j++)
			{ 
				d=fabs(iR1[i][j]);
				if ((i != j)&&(d > fm))
				{ 
					fm = d; 
					p = i;
					q = j;
				}
			}
		}

        if (fm < eps)
		{
			for (i = 0; i < iNum; ++i)
				EigenValue[i] = iR1[i][i];
			l = 200;
		}
        
		l=l + 1;
        x=-iR1[p][q];
		y=(iR1[q][q]-iR1[p][p]) / 2.0;
        omega=x / sqrt(x * x + y * y);

        if (y<0.0) 
			omega=-omega;

        sn=1.0+sqrt(1.0-omega*omega);
        sn=omega / sqrt(2.0 * sn);
        cn=sqrt(1.0-sn*sn);
        fm=iR1[p][p];
        iR1[p][p]=fm*cn*cn+iR1[q][q]*sn*sn+iR1[p][q]*omega;
        iR1[q][q]=fm*sn*sn+iR1[q][q]*cn*cn-iR1[p][q]*omega;
        iR1[p][q]=0.0;
		iR1[q][p]=0.0;
        for (j=0; j<=iNum-1; j++)
		{
			if ((j!=p)&&(j!=q))
			{
				fm=iR1[p][j];
				iR1[p][j]=fm*cn+iR1[q][j]*sn;
				iR1[q][j]=-fm*sn+iR1[q][j]*cn;
			}
		}

        for (i=0; i<=iNum-1; i++)
		{
			if ((i!=p)&&(i!=q))
            { 
				fm=iR1[i][p];
				iR1[i][p]=fm*cn+iR1[i][q]*sn;
				iR1[i][q]=-fm*sn+iR1[i][q]*cn;
            }
		}

        for (i = 0; i <= iNum-1; i++)
        { 
            fm=EigenVector[i][p];
            EigenVector[i][p]=fm*cn+EigenVector[i][q]*sn;
            EigenVector[i][q]=-fm*sn+EigenVector[i][q]*cn;
        }
    }
    
	printf("特征值：\n");
	for (i = 0; i < iNum; ++i)
	{
		EigenValue[i] = iR1[i][i];
        printf("%f ",EigenValue[i]);
	}
	printf("\n特征向量：\n");

	for(i = 0;i < iNum;i++)
	{
		for(j = 0;j < iNum;j++)
			printf("%f ",EigenVector[i][j]);
		printf("\n");
 	}
 	return true;
}