特征值问题数值解

介绍

当使用FEM求解特征值问题时,一般可以表述为:

在大多数工程问题中,A,B一般是对称矩阵;如果要求的问题是结构自由振动,那么A对应刚度矩阵K;B对应质量矩阵M; \lambda 是特征值; \vec{X} 是特征向量或称为模态振型.假设A,B已经引入边界条件.

上式可重写为:

有解:

可知,当且仅当系数矩阵行列式为0,方程才有非0解,即:

将行列式展开可得到关于 \lambda 的n次多项式,其根就是特征值 \lambda_i ,把 \lambda_i 带入方程可以求得 \vec{X}_i .这称为特征对.根据矩阵论, X_j 是特征向量,那么n个特征向量的线性组合也是特征向量.

在商软中,有一种归一化方法是,每个特征向量的最大值归一化为1.0,这称为位移归一化.如下:

另一种归一化称为质量归一,满足:

标准特征值问题

特征值问题分为广义特征值问题和标准特征值问题.广义特征值问题形如:

标准特征值问题形如:

求解广义特征值的第一步是转换为标准特征值问题:

[A]X=\lambda [B]X \Rightarrow [B]^{-1}[A]X=\lambda X \tag{1}

虽然A,B一般是对称矩阵,但 H=B^{-1}A 一般则不对称.可以利用choleski分解得到对称的H矩阵.假设B是对称正定矩阵,则 B=U^TU

B带入公式1后:

定义 \vec{Y}=U\vec{X} ,则有:

([H]-\lambda[I]) \vec{Y}=\overrightarrow{0}\tag{2}

[H]是对称矩阵,即:

[H]=\left([U]^{T}\right)^{-1}[A][U]^{-1}

求解公式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

计算最大特征值

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 ,【这可以作为停止迭代的标志】,进一步:

\lim_{p\rightarrow \inf}(\frac{[H]^p\vec{Z_0}}{\lambda_1^p})=a_1\vec{X_1}

\lambda_1 也可以通过瑞利商R得到:

如果 [H]\vec{X}=\lambda\vec{X} ,R则为 \lambda .进一步,第i次迭代中计算R:

只要相邻两次迭代的R相差小于tol 值,那么 \lambda_1=R_i

计算最小特征值

我们可以通过计算下式来获得最小特征值与相应的特征向量:

计算方法: 1)左乘H的逆,可得:

2)改写为:

这意味着可以像计算最大特征值那样计算[H]的最小特征值.

注意:计算最小特征值之前,要先找到[H]的逆,虽然会有额外的计算,但在某些情况下是比较好的办法.

计算中间特征值

对 \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 同堂可以得到.按照类似的过程循环接下去就可以得到其他特征值.

瑞利里兹-子空间迭代法

这是一个求解大规模自由度系统的少量特征值的有效算法;下面给出算法的简要介绍,详细细节见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.

Step01

取q个初始迭代向量如下,其中p是要求的特征值个数. Bathe and Wilson发现q=min(2p,p+8)可以获得较好的收敛性;

定义初始模态矩阵X0:

此时迭代次数k=0.

T.C. Cheu, C.P. Johnson, R.R. Craig Jr., Computer algorithms for calculating efﬁcient initial vectors for subspace iteration method, International Journal for Numerical Methods in Engineering 24 (1987) 1841-1848.给出了为子空间迭代法计算有效的初始向量的一种方法

Step02

现在需要利用子空间迭代法生成改进模态矩阵 X_{k+1} :

X_{k+1} 可由下式计算:

2)进一步计算:

3)计算缩减系统的特征值和特征向量,并求得 Q_{k+1} 和 \lambda_{k+1} :

[A_{k+1}][Q_{k+1}]=[B_{k+1}][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个最小特征值预测的近似值上测量.

返回主页

目录