目的
求一個實對稱矩陣的所有特征值和特征向量。
前置知識
對於一個實對稱矩陣\(A\)，必存在對角陣\(D\)和正交陣\(U\)滿足$$D=U^TAU$$\(D\)的對角線元素為\(A\)的特征值，\(U\)的列向量為\(A\)的特征向量。
定義\(n\)階旋轉矩陣$$G(p,q,\theta)=
\begin{bmatrix}
1 & & & & & \cdots& & & & & 0\
&\ddots & & & & & & & & &\
& &1 & & & & & & & &\
& & &\cos\theta & & & &-\sin\theta & & &\
& & & &1 & &0 & & & &\
& & & & &\ddots & & & & &\
& & & &0 & &1 & & & &\
& & &\sin\theta & & & &\cos\theta & & &\
& & & & & & & &1 & &\
& & & & & & & & &\ddots &\
0& & & & & & & &0 & &1
\end{bmatrix}
\[即在單位矩陣的基礎上，修改$a_{pp}=a_{qq}=\cos\theta,a_{qp}=-a_{pq}=\sin\theta$ 對於$n$階向量$\alpha$，$\alpha\cdot G(p,q,\theta)$的幾何意義是把$\alpha$在與第$p$維坐標軸和第$q$維坐標軸平行的平面內旋轉角度$\theta$，並且旋轉后的模長保持不變。 ## 算法原理 大概思路使通過旋轉變換使非對角線上的元素不斷變小，最后得到與原矩陣相似的對角矩陣。 每次找到矩陣$A$絕對值最大的的非對角線元素，設為$a_{pq}$，令$U=G(p,q,\theta)$，將$A$變換為$U^TAU$ 變換后的值為  通過令$b_{p,q}=0$解得$$\theta=\frac{1}{2}\arctan\frac{2a_{pq}}{a_{qq}-a_{pp}}$$特別地當$a_{qq}=a_{pp}$時$\theta=\frac{\pi}{4}$ 注意到旋轉操作並不會改變每個行向量或列向量的模長，即矩陣$A$的F-范數$||A||_F=\sqrt{\sum_i\sum_ja_{ij}^2}$是不變的，並且通過計算可以得出$$b_{ip}^2+b_{iq}^2=a_{ip}^2+a_{iq}^2$$從而可以得知非對角線元素的平方和變小，對角線上元素的平方和增大，故非主對角線上元素的平方和收斂。 ## 算法流程 （1）令矩陣$T=E$，即初始化單位矩陣 （2）找到$A$中絕對值最大的非對角選元素$a_{pq}$ （3）找到對應的角度$\theta$，構造矩陣$U=G(p,q,\theta)$ （4）令$A=U^TAU,T=TU$ （5）不停地重復（2）到（4），直到$a_{pq}<\epsilon$或迭代次數超過某個限定值，則$A$的對角線元素近似等於$A$的特征值，$T$的列向量為$A$的特征向量 ## 代碼 ```c #include<bits/stdc++.h> using namespace std; const int N=1005; const double eps=1e-5; const int lim=100; int n,id[N]; double key[N],mat[N][N],EigVal[N],EigVec[N][N],tmpEigVec[N][N]; bool cmpEigVal(int x,int y) { return key[x]>key[y]; } void Find_Eigen(int n,double (*a)[N],double *EigVal,double (*EigVec)[N]) { for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) EigVec[i][j]=0; for (int i=1;i<=n;i++) EigVec[i][i]=1.0; int count=0; while (1) { //統計迭代次數 count++; //找絕對值最大的元素 double mx_val=0; int row_id,col_id; for (int i=1;i<n;i++) for (int j=i+1;j<=n;j++) if (fabs(a[i][j])>mx_val) mx_val=fabs(a[i][j]),row_id=i,col_id=j; if (mx_val<eps||count>lim) break; //進行旋轉變換 int p=row_id,q=col_id; double Apq=a[p][q],App=a[p][p],Aqq=a[q][q]; double theta=0.5*atan2(-2.0*Apq,Aqq-App); double sint=sin(theta),cost=cos(theta); double sin2t=sin(2.0*theta),cos2t=cos(2.0*theta); a[p][p]=App*cost*cost+Aqq*sint*sint+2.0*Apq*cost*sint; a[q][q]=App*sint*sint+Aqq*cost*cost-2.0*Apq*cost*sint; a[p][q]=a[q][p]=0.5*(Aqq-App)*sin2t+Apq*cos2t; for (int i=1;i<=n;i++) if (i!=p&&i!=q) { double u=a[p][i],v=a[q][i]; a[p][i]=u*cost+v*sint;a[q][i]=v*cost-u*sint; u=a[i][p],v=a[i][q]; a[i][p]=u*cost+v*sint;a[i][q]=v*cost-u*sint; } //計算特征向量 for (int i=1;i<=n;i++) { double u=EigVec[i][p],v=EigVec[i][q]; EigVec[i][p]=u*cost+v*sint;EigVec[i][q]=v*cost-u*sint; } } //對特征值排序 for (int i=1;i<=n;i++) id[i]=i,key[i]=a[i][i]; std::sort(id+1,id+n+1,cmpEigVal); for (int i=1;i<=n;i++) { EigVal[i]=a[id[i]][id[i]]; for (int j=1;j<=n;j++) tmpEigVec[j][i]=EigVec[j][id[i]]; } for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) EigVec[i][j]=tmpEigVec[i][j]; //特征向量為列向量 } int main() { scanf("%d",&n); for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) scanf("%lf",&mat[i][j]); Find_Eigen(n,mat,EigVal,EigVec); printf("EigenValues = "); for (int i=1;i<=n;i++) printf("%lf ",EigVal[i]); printf("\nEigenVector =\n"); for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) printf("%lf%c",EigVec[i][j],j==n?'\n':' '); return 0; } ```\]