一 算法原理

雅可比方法用于求解实对称矩阵的特征值和特征向量,对于实对称矩阵$A$,必有正交矩阵$U$,使得$U^{T}AU=D$.$D$是一个对角阵,主对角线的元素是矩阵$A$的特征值,正交矩阵$U$的每一列对应于属于矩阵$D$的主对角线对应元素的特征向量.
雅可比方法用平面旋转对矩阵$A$做相似变换,化$A$为对角阵,从而求解出特征值和特征向量.旋转矩阵$U_p{}_q$,是一个单位阵在第$p$行,第$p$列,第$q$行,第$q$列,元素为$cos\phi$,第$p$行第$q$列为$-sin\phi$,第$q$行第$p$列为$sin\phi$.对于这样的平面旋转矩阵,不难验证其是一种正交矩阵.因此对于向量$x$,$U_p{}_{q}x$等同于把第$p$个坐标轴和第$q$个坐标轴共同所确定的平面旋转了$\phi$度.记矩阵$A_1=U_p{{}_q}^{T}A{U}_p{}_q$.因为旋转矩阵是正交阵,因此实际上矩阵$A_1$与矩阵$A$是相似的,因此其特征值是相同的.
设矩阵$A_1$第$i$行,第$j$列的元素为$b_i{}_j$,矩阵$A$的第$i$行,第$j$列的元素为$a_i{}_j$($i=0,1,2,...,n-1,j=0,1,2,...,n-1$).式(1-1-1)给出了两矩阵元素之间的运算关系.
$$\{\begin{array}{l}{b}_p{}_p=a_p{}_{p}co{s}^2\phi +a_q{}_{q}si{n}^2\phi +2a_p{}_{q}cos\phi sin\phi \\ b_q{}_q=a_p{}_{p}si{n}^2\phi +a_q{}_{q}co{s}^2\phi -2a_p{}_{q}cos\phi sin\phi \\ b_p{}_q=b_q{}_p=\frac{1}{2}(a_q{}_q-a_p{}_p)sin2\phi +a_p{}_{q}cos2\phi \\ b_p{}_i=a_p{}_{i}cos\phi +a_q{}_{i}sin\phi ,(i\ne p,q)\\ b_q{}_i=-a_p{}_{i}sin\phi +a_q{}_{i}cos\phi ,(i\ne p,q)\\ b_j{}_p=a_j{}_{p}cos\phi +a_j{}_{q}sin\phi ,(j\ne p,q)\\ b_j{}_q=-a_j{}_{q}sin\phi +a_j{}_{q}cos\phi ,(j\ne p,q)\\ b_i{}_j=b_j{}_i=a_i{}_j,i\ne p,q;j\ne p,q\end{array}\text{\hspace{1em}(1-1-1)}$$
其中有两点需要说明:(1) $p$,$q$分别是前一次的迭代矩阵的非主对角线上绝对值最大元素的行列号
(2)$\phi$是旋转角度,可以由式(1-1-2)确定
$$tan2\phi =\frac{-2a_p{}_q}{a_q{}_q-a_p{}_p}\text{\hspace{1em}(1-1-2)}$$
归纳得到雅可比方法求解矩阵特征值和特征向量的具体步骤如下:
(1).初始化特征向量为对角阵V,主对角线元素为1,其他元素为0.
(2).在$A$的非主对角线的元素中,找到绝对值最大元素$a_p{}_q$.
(3).用式(1-1-2)计算出旋转矩阵
(4).计算矩阵$A1$,用当前的矩阵$V$乘旋转矩阵得到当前的特征矩阵$V$.
(5).若当前迭代的矩阵$A$的非主对角线元素最大值小于给定阈值,停止计算,否则执行上述过程.停止计算时,特征值为矩阵$A$的主对角线元素,特征矩阵为矩阵$V$.

二 C语言实现

#include<stdio.h>
#include<stdlib.h>
float** Matrix_Jac_Eig(float **array, int n, float *eig);
int Matrix_Free(float **tmp, int m, int n);
int main(void)
{
	int n;
	printf("请输入矩阵维度:\n");
	scanf("%d", &n);
	float **array = (float **)malloc(n * sizeof(float *));
	if (array == NULL)
	{
		printf("error :申请数组内存空间失败\n");
		return -1;
	}
	for (int i = 0; i < n; i++)
	{
		array[i] = (float *)malloc(n * sizeof(float));
		if (array[i] == NULL)
		{
			printf("error :申请数组子内存空间失败\n");
			return -1;
		}
	}
	printf("请输入矩阵元素:\n");
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			scanf("%f", &array[i][j]);
		}
	}
	float *eig = (float *)malloc(n * sizeof(float));
	float **Result = Matrix_Jac_Eig(array, n, eig);
	printf("特征矩阵元素:\n");
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			printf("%f ", Result[i][j]);
		}
		printf("\n");
	}
	printf("特征矩阵元素:\n");
	for (int i = 0; i < n; i++)
	{
		printf("%f \n", eig[i]);
	}
	Matrix_Free(Result, n, n);
	free(eig);
	eig = NULL;
	return 0;
}
float** Matrix_Jac_Eig(float **array, int n, float *eig)
{
	//先copy一份array在temp_mat中，因为我实在堆区申请的空间,在对其进行处理
	//的过程中会修改原矩阵的值,因此要存储起来,到最后函数返回的
	//时候再重新赋值
	int i, j, flag, k;
	flag = 0;
	k = 0;
	float sum = 0;
	float **temp_mat = (float **)malloc(n * sizeof(float *));
	for (i = 0; i < n; i++)
	{
		temp_mat[i] = (float *)malloc(n * sizeof(float));
	}
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			temp_mat[i][j] = array[i][j];
		}
	}
	//判断是否为对称矩阵
	for (i = 0; i < n; i++)
	{
		for (j = i; j < n; j++)
		{
			if (array[i][j] != array[j][i])
			{
				flag = 1;
				break;
			}
		}
	}
	if (flag == 1)
	{
		printf("error in Matrix_Eig: 输入并非是对称矩阵:\n");
		return NULL;
	}
	else
	{
		//开始执行算法
		int p, q;
		float thresh = 0.0000000001;
		float max = array[0][1];
		float tan_angle, sin_angle, cos_angle;
		float **result = (float **)malloc(n * sizeof(float *));
		if (result == NULL)
		{
			printf("error in Matrix_Eig:申请空间失败\n");
			return NULL;
		}
		float **result_temp = (float **)malloc(n * sizeof(float *));
		if (result_temp == NULL)
		{
			printf("error in Matrix_Eig:申请空间失败\n");
			return NULL;
		}
		float **rot = (float **)malloc(n * sizeof(float *));
		if (rot == NULL)
		{
			printf("error in Matrix_Eig:申请空间失败\n");
			return NULL;
		}
		float **mat = (float **)malloc(n * sizeof(float *));
		if (mat == NULL)
		{
			printf("error in Matrix_Eig:申请空间失败\n");
			return NULL;
		}
		for (i = 0; i < n; i++)
		{
			result[i] = (float *)malloc(n * sizeof(float));
			if (result[i] == NULL)
			{
				printf("error in Matrix_Eig:申请子空间失败\n");
				return NULL;
			}
			result_temp[i] = (float *)malloc(n * sizeof(float));
			if (result_temp[i] == NULL)
			{
				printf("error in Matrix_Eig:申请子空间失败\n");
				return NULL;
			}
			rot[i] = (float *)malloc(n * sizeof(float));
			if (rot[i] == NULL)
			{
				printf("error in Matrix_Eig:申请子空间失败\n");
				return NULL;
			}
			mat[i] = (float *)malloc(n * sizeof(float));
			if (mat[i] == NULL)
			{
				printf("error in Matrix_Eig:申请子空间失败\n");
				return NULL;
			}
		}
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				if (i == j)
				{
					result[i][j] = 1;
				}
				else
				{
					result[i][j] = 0;
				}
			}
		}
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				if (i == j)
				{
					mat[i][j] = 1;
				}
				else
				{
					mat[i][j] = 0;
				}
			}
		}
		max = array[0][1];
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				if (i == j)
				{
					continue;
				}
				else
				{
					if (fabs(array[i][j]) >= fabs(max))
					{
						max = array[i][j];
						p = i;
						q = j;
					}
					else
					{
						continue;
					}
				}
			}
		}
		while (fabs(max) > thresh)
		{
			if (fabs(max) < thresh)
			{
				break;
			}
			tan_angle = -2 * array[p][q] / (array[q][q] - array[p][p]);
			sin_angle = sin(0.5*atan(tan_angle));
			cos_angle = cos(0.5*atan(tan_angle));
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{

					if (i == j)
					{
						mat[i][j] = 1;
					}
					else
					{
						mat[i][j] = 0;
					}
				}
			}
			mat[p][p] = cos_angle;
			mat[q][q] = cos_angle;
			mat[q][p] = sin_angle;
			mat[p][q] = -sin_angle;
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					rot[i][j] = array[i][j];
				}
			}
			for (j = 0; j < n; j++)
			{
				rot[p][j] = cos_angle*array[p][j] + sin_angle*array[q][j];
				rot[q][j] = -sin_angle*array[p][j] + cos_angle*array[q][j];
				rot[j][p] = cos_angle*array[j][p] + sin_angle*array[j][q];
				rot[j][q] = -sin_angle*array[j][p] + cos_angle*array[j][q];
			}
			rot[p][p] = array[p][p] * cos_angle*cos_angle +
				array[q][q] * sin_angle*sin_angle +
				2 * array[p][q] * cos_angle*sin_angle;
			rot[q][q] = array[q][q] * cos_angle*cos_angle +
				array[p][p] * sin_angle*sin_angle -
				2 * array[p][q] * cos_angle*sin_angle;
			rot[p][q] = 0.5*(array[q][q] - array[p][p]) * 2 * sin_angle*cos_angle +
				array[p][q] * (2 * cos_angle*cos_angle - 1);
			rot[q][p] = 0.5*(array[q][q] - array[p][p]) * 2 * sin_angle*cos_angle +
				array[p][q] * (2 * cos_angle*cos_angle - 1);
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					array[i][j] = rot[i][j];
				}
			}
			max = array[0][1];
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					if (i == j)
					{
						continue;
					}
					else
					{
						if (fabs(array[i][j]) >= fabs(max))
						{
							max = array[i][j];
							p = i;
							q = j;
						}
						else
						{
							continue;
						}
					}
				}
			}
			for (i = 0; i < n; i++)
			{
				eig[i] = array[i][i];
			}
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					sum = 0;
					for (k = 0; k < n; k++)
					{
						sum = sum + result[i][k] * mat[k][j];
					}
					result_temp[i][j] = sum;
				}
			}
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					result[i][j] = result_temp[i][j];
				}
			}
		}
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
			{
				array[i][j] = temp_mat[i][j];
			}
		}
		Matrix_Free(result_temp, n, n);
		Matrix_Free(rot, n, n);
		Matrix_Free(mat, n, n);
		Matrix_Free(temp_mat, n, n);
		return result;
	}
}
int Matrix_Free(float **tmp, int m, int n)
{
	int i, j;
	if (tmp == NULL)
	{
		return(1);
	}
	for (i = 0; i < m; i++)
	{
		if (tmp[i] != NULL)
		{
			free(tmp[i]);
			tmp[i] = NULL;
		}
	}
	if (tmp != NULL)
	{
		free(tmp);
		tmp = NULL;
	}
	return(0);
}

三 结果

1. 2阶实对称矩阵特性

定理：2阶实对称矩阵H的特征值是实数

H=[a,b;b,c]    

    a,b,c是实数，λ 是特征值
A=[a-λ,b;b,c-λ]    

    特征值求解方法为：(a- λ )(c- λ) - b2  = 0

    求解方程得到两个根为：
λ=（a+c）±（a+c）2-4（ac-b2）2  

                   （a+c）2-4ac-b2=a-c2+4b2≥0

    所以，在a、b、c为实数时，特征值也是实数。

2、特征向量

根据特征值和特征向量的定义：HX=λX，（H-λE）X = 0；因此方程若有解，则

det（H-λE）=0；

设
A=[a-λ,b;b,c-λ]    

则有-b/(a-λ) = (c-λ)/b,   线性齐次方程组AX=0有非零解，其中之一解向量 [1,-b/(a-λ)],归一化后得到标准解。