Note_NAA_特征值问题数值解
介绍
当使用FEM求解特征值问题时，一般可以表述为：
在大多数工程问题中，A,B一般是对称矩阵；如果要求的问题是结构自由振动，那么A对应刚度矩阵K；B对应质量矩阵M；\(\lambda\)是特征值；\(\vec{X}\)是特征向量或称为模态振型。假设A,B已经引入边界条件，参见Note_Fem边界条件的处理和实现。
上式可重写为：
有解：
可知，当且仅当系数矩阵行列式为0，方程才有非0解，即：
将行列式展开可得到关于\(\lambda\)的n次多项式，其根就是特征值\(\lambda_i\),把\(\lambda_i\)带入方程可以求得\(\vec{X}_i\).这称为特征对。根据矩阵论，\(X_j\)是特征向量，那么n个特征向量的线性组合也是特征向量。
在商软中，有一种归一化方法是，每个特征向量的最大值归一化为1.0，这称为位移归一化。如下：
另一种归一化称为质量归一，满足：
标准特征值问题
特征值问题分为广义特征值问题和标准特征值问题。广义特征值问题形如：
标准特征值问题形如：
求解广义特征值的第一步是转换为标准特征值问题：
虽然A,B一般是对称矩阵，但\(H=B^{-1}A\)一般则不对称。可以利用choleski分解得到对称的H矩阵。假设B是对称正定矩阵，则\(B=U^TU\)
B带入公式1后：
定义\(\vec{Y}=U\vec{X}\),则有：
[H]是对称矩阵，即：
求解公式2中的\(\lambda_i\)和\(\vec{Y_i}\)后，\(\vec{X_i}=U^{-1}\vec{Y}\).
标准特征值的数值解
求解方法分类：
- transformation methods
- Jacobi method
- Givens method
- Householder method
- etc
- iterative methods
- power method
- etc
transformation methods适合用于求所有特征值对；迭代法适合求解少数几个特征值对。
Jacobi Method
solving:
Jacobi方法基于： 实对称矩阵H只有实特征值\(\lambda_i\)，且存在实正交矩阵P，P满足\(P^THP=diag\{\lambda_i\}\),P矩阵的列是特征向量。
P矩阵可以由一系列形如下式的旋转矩阵的乘积而得到。
\([P_1]\)中除了第i列和第j行中的元素外，所有元素都与单位矩阵\(I\)中的元素相同。如果Sin值和Cos值出现在(i, i)，(i, j), (j, i),(j, j), 那\([P_1]^TH[P]\)的对应元素为：
\(\theta\)满足：
然后令\(h_{ij}=h_{ji}=0\).因此在Jacobi法每一步中都会令一对非对角元素为0.0；在下一步中，虽然该方法减少了一对新的零，但它引入了非零的贡献到以前的零位置。在下一步中使用下入下式的矩阵[H_i]代替[H]：
关于每一步的i,j选择：每一步开始时，可以选取[H_i]非对角元素的最大值所在的i，j；一直迭代到H_i矩阵非对角元素全为0.0，或小于某个阈值。
最后的特征向量矩阵\(P=P_1P_2P_3...\)
以上时jacobi法的基本操作，由此衍生了循环jacobi，过关jacobi等方法。
https://www.cnblogs.com/tlnshuju/p/6726000.html https://zhuanlan.zhihu.com/p/262879394 https://blog.csdn.net/qq_43133135/article/details/126502543
Power Method
1. 计算最大特征值
power method是求解最大特征值或主特征值(\(\lambda_1\))和对应特征向量\(\vec{X_1}\)的最简单的迭代法。假设H为实对称矩阵，存在N个独立特征向量；
首先，选择初始向量\(\vec{Z_0}\),通过迭代生成一系列向量\(\vec{Z_1},\vec{Z_2}....\):
因此，p-th向量有：
在第p次迭代会得到\(Z_p\),如果满足：
则停止迭代。\(z_{p,j}, z_{p-1,j}\)表示\(Z_p\)和\(Z_{p-1}\)的第j个分量。\(\lambda_1\)是要求的特征值；
收敛性解释如下
Z0向量可以表示为特征向量得到线性组合：
如果\(\lambda_i\)是\(X_i\)对应的特征值，则：
根据以上，当p=+inf,则：
因此取\(H^pZ_0\)和\(H^{p-1}Z_0\)的任一分量相除，应该有相同的极限值\(\lambda_1\)，【这可以作为停止迭代的标志】，进一步：
\(\lambda_1\)也可以通过瑞利商R得到：
如果\([H]\vec{X}=\lambda\vec{X}\),R则为\(\lambda\).进一步，第i次迭代中计算R：
只要相邻两次迭代的R相差小于tol 值，那么\(\lambda_1=R_i\)
2. 计算最小特征值
我们可以通过计算下式来获得最小特征值与相应的特征向量：
计算方法： 1）左乘H的逆，可得：
2）改写为：
这意味着可以像计算最大特征值那样计算[H]的最小特征值。
注意：计算最小特征值之前，要先找到[H]的逆，虽然会有额外的计算，但在某些情况下是比较好的办法。
3. 计算中间特征值
对\(\vec{X_1}\)进行归一化，令X1的第一个分量\(x_1=1\),即：
记\(\vec{r}^T\)为[H]的第一行元素，\(\vec{r}^T=\{h_{11},h_{12},...,h_{1n},\}\),那么有：
令另一个主特征值为\(\lambda_2\),归一化\(\vec{X_2}\),令\(\vec{X_2}\)的第一个分量为1.0：
eq7.56表明：\(\lambda_2\)是\([H]-\left [ \underset{\sim }{H} \right ]\)的特征值；\([H]-\left [ \underset{\sim }{H} \right ]\)的第一行全为0，而X2-X1的第一个元素为0；通过删去\([H]-\left [ \underset{\sim }{H} \right ]\)的第一行和第一列可以获得\([H_2]\),然后，计算\([H_2]\)的主特征值和特征向量，并且把特征向量的第一个元素置0，可以获得Z1.可以发现X2-X1必然是Z1的倍数,即：
乘子a为：
至此，特征向量X2就可以得到了，\(\lambda_2\)同堂可以得到。按照类似的过程循环接下去就可以得到其他特征值。
瑞利里兹-子空间迭代法
这是一个求解大规模自由度系统的少量特征值的有效算法；下面给出算法的简要介绍，详细细节见1
Step01
取q个初始迭代向量如下，其中p是要求的特征值个数。Bathe and Wilson发现q=min(2p,p+8)可以获得较好的收敛性；
定义初始模态矩阵X0：
此时迭代次数k=0.
Cheu2给出了为子空间迭代法计算有效的初始向量的一种方法
Step02
现在需要利用子空间迭代法生成改进模态矩阵\(X_{k+1}\)： 1）\(X_{k+1}\)可由下式计算：
2）进一步计算：
3）计算缩减系统的特征值和特征向量，并求得\(Q_{k+1}\)和\(\lambda_{k+1}\)：
4）找到原始系统特征向量的改进近似值：
It is assumed that the iteration vectors converging to the exact eigenvectors,\(\vec{X_1^(extra)},\vec{X_2^(extra)}.....\) ,并且就是[X_{k+1}]的列向量 It is assumed that the vectors in [X0] are not orthogonal to one of the required eigenvectors
Step03
如果相邻两次迭代过程k-1和k中求得的近似特征值\(\lambda_i^{(k)} and\lambda_i^{(k+1)}\)满足下式则认为收敛：
\(\epsilon\)=1e-6;虽然迭代是用q个向量（q > p）执行的，但收敛性仅在p个最小特征值预测的近似值上测量.
-
K.J. Bathe, E.L. Wilson, Large eigenvalue problems in dynamic analysis, Journal of Engineering Mechanics Division, Proceedings of ASCE 98 (EM6) (1972) 1471-1485. ↩
-
T.C. Cheu, C.P. Johnson, R.R. Craig Jr., Computer algorithms for calculating efficient initial vectors for subspace iteration method, International Journal for Numerical Methods in Engineering 24 (1987) 1841-1848. ↩