多自由度振动
约 6034 个字 3 张图片 预计阅读时间 40 分钟
Notes
多自由度系统一定会出一道大题，如果是系统计算分析的题目一般是两个自由度或三个自由度，多的不好算
双自由度系统
我们之前所关注的都是单自由度系统的振动问题，但是实际生活中有很多更为复杂的振动，不能单纯化简为单自由度系统，必须化成双自由度甚至多自由度系统的问题来解决。随着振系自由度数目的增多，振动问题求解的工作量越来越繁重。不过，多自由度振系的许多基本概念都可以通过双自由度振系的问题来说明；而且，关于双自由度振系的理论本身有很重要。
考虑如下双自由度振动系统：
考虑一种简单地特殊情况：取 \(k_3 = k_1,c_3=c_1\) 则可以写出系统运动方程
\[
\begin{aligned}
&m \ddot{x}_1= F_1(t) + k_2(x_2-x_1)-k_1x_1+c_2(\dot{x}_2-\dot{x}_1)-c_1\dot{x}_1\\
&m \ddot{x}_2= F_2(t) - k_2(x_2-x_1)-k_3x_2-c_2(\dot{x}_2-\dot{x}_1)-c_2\dot{x}_1
\end{aligned}
\]
移项并写成矩阵形式有
\[
\begin{bmatrix}
m&0\\
0&m
\end{bmatrix}
\begin{bmatrix}
\ddot{x}_1\\
\ddot{x}_2
\end{bmatrix}+
\begin{bmatrix}
c_1+c_2&-c_2\\
-c_2&c_1+c_2
\end{bmatrix}
\begin{bmatrix}
\dot{x}_1\\
\dot{x}_2
\end{bmatrix}+
\begin{bmatrix}
k_1+k_2&-k_2\\
-k_2&k_1+k_2
\end{bmatrix}
\begin{bmatrix}
{x}_1\\
{x}_2
\end{bmatrix}=
\begin{bmatrix}
F_1(t)\\
F_2(t)
\end{bmatrix}
\]
记作
\[
M\ddot{X}+C\dot{X}+KX=F
\]
先考虑保守系统，此时
\[
M\ddot{X} +KX= 0
\]
假定存在
\[
X= X_0 \sin(\omega t+\varphi)
\]
带入消去 \(\sin(\omega t +\varphi)\)
\[
(K-\omega^2M)X=0
\]
为了满足方程有解，此时行列式 \(\det[K-\omega^2M]\) 应为零，对于我们所研究的系统，有
\[
(k_1+k_2-\omega^2m)^2=k_2^2
\]
解得
\[
\omega_1^2= \frac{k_1}{m}\quad w^2_2=\frac{k_1+2k_2}{m}
\]
对于 \(\omega_1\) 带入可解得振幅满足关系 \(X_{10}= X_{20}\)，定义振型，事实上它就是特征值对应的特征向量
\[
\varphi_1= \begin{bmatrix}
1\\
1
\end{bmatrix}\quad(\text{归一化形式：}\begin{bmatrix}
\frac{1}{\sqrt{2}}\\
\frac{1}{\sqrt{2}}
\end{bmatrix})
\]
对于 \(\omega_2\) 带入可解得振幅满足关系 \(X_{10}= -X_{20}\)，振型
\[
\varphi_2= \begin{bmatrix}
1\\
-1
\end{bmatrix}\quad(\text{归一化形式：}\begin{bmatrix}
\frac{1}{\sqrt{2}}\\
-\frac{1}{\sqrt{2}}
\end{bmatrix})
\]
显然 \(\varphi_1\) 与 \(\varphi_2\) 正交
在振动过程中，每个质量振动所对应的实际幅值是不确定的，但是各个质量的振动幅度与特征向量中的各个元素成比例。振动系统的固有频率取决于质量矩阵和刚度矩阵，主振型也是由质量矩阵和刚度矩阵决定的，是系统固有的振动性质。
线性代数告诉我们，不同特征值所对应的特征向量之间相互正交，也就是说，多自由振动系统的主振型相互正交，利用其正交性我们有
\[
\varphi_i^TM\varphi_j = \left\{\begin{align}
&0&i\neq j\\
&\overline{M}_i(2m) &i=j
\end{align}\right.
\]
\[
\varphi_i^TK\varphi_j = \left\{\begin{align}
&0&i\neq j\\
&\overline{K}_1 = 2k_1 &i=j=1\\
&\overline{K}_2 = 2k_1+4k_2 &i=j=2
\end{align}\right.
\]
所以
$$
\omega_i^2= \frac{\overline{K}_i}{\overline{M}_i}
$$
我们给出振幅曲线和对应简化的物理模型：
- 对于 \(\varphi_1\) 我们可以视作两个质量以相同的振幅同向运动，中间弹簧无变形。
- 对于 \(\varphi_2\) 我们可以视作两个质量以相同的振幅反向向运动，中间弹簧中点不动。
对于其他情况（两端弹簧质量不相同），我们仍可以用相同的方法求得两个质量的振幅比 \(\frac{X_{10}}{X_{20}}\) 关系。此时仍会存在一正一负两种振幅比 \(\mu_1,\mu_2\) ，称\(\mu_1>0\)对应的振动为第一振型（或低振型振动）称\(\mu_2<0\)对应的振动为第二振型（或高振型振动），如下图所示：
回到之前的方程，我们尝试对耦合的运动方程解耦：
取
$$
X=\begin{bmatrix}
\varphi_1&\varphi_2
\end{bmatrix}\begin{bmatrix}
q_1(t)\
q_2(t)
\end{bmatrix}=\begin{bmatrix}
1\1
\end{bmatrix}q_1(t)+\begin{bmatrix}
1\-1
\end{bmatrix}q_2(t)
$$
带入后仍是个矩阵形式，我们通过左乘 \(\varphi_1^T\)或 \(\varphi^T_2\)化成标量方程可得
\[
\begin{align}
&\varphi^T_1:2m \ddot{q}_1+2c_1 \dot{q}+2k_1q_1 = F_1+F_2\\
&\varphi^T_2:2m \ddot{q}_2+(2c_1+4c_2)\dot{q}+(2k_1+4k_2)\dot{q}_2 = F_1-F_2
\end{align}
\]
但此时要求阻尼为比例阻尼，即
\[
C=\alpha M +\beta K
\]
不然难以解耦。
多自由振动系统
对于多自由振动系统，直接利用牛顿力学建立运动方程是一件非常复杂的事情。我们考虑利用Lagrange方程建立系统方程。
给出 Lagrange 作用量
\[
L=T-V = \frac{1}{2}m\dot{q}_i\dot{q}_j - \frac{1}{2}R_{ij}q_jq_i
\]
带入 Lagrange 方程中有：
\[
\frac{d}{dt}\left[\frac{\partial L}{ \partial \dot{q}_i}\right]-\frac{ \partial L}{\partial q_i} =Q_i(t)
\]
\(Q_i(t)\) 为广义力，它包括阻尼力和外加激励。
如果在 Lagrange 量中引入能量耗散函数
\[
D =\frac{1}{2}c_{ij}\dot{q}_i\dot{q}_j
\]
则将 Lagrange 方程改写成
\[
\frac{d}{dt}\left[\frac{\partial L}{ \partial \dot{q}_i}\right]-\frac{ \partial L}{\partial q_i} + \frac{\partial D }{\partial q_i}=Q_i(t)
\]
此时广义力只有外加激励的作用结果。
将上述的 Lagrange 方程展开可得多自由度振动系统的振动微分方程
\[
M\ddot{q}+C\dot{q}+Kq=F
\]
这个方程对于单自由振动系统一定成立，但是对于多自由度振动系统，系统中仍存在其他作用因素
- 科氏力（也称陀螺力） \(G\dot{q}\)
- 因张力刚度发生的变化 \(\overset{\sim}{K} q\)
这些因素在后面暂时不会考虑，我们对科氏力做一个简要说明，首先系数矩阵 \(G\) 是一个反对称矩阵
\[
G= -G^T
\]
所以无法将其与阻力系数矩阵 \(C\) 合并。
科氏力是一个保守力，考虑一个陀螺，陀螺转动能保持长时间的稳定，我们可以将其视作是一种保守系统，即没有能量的耗散。
保守状态的多自由振动系统
模仿双自由度系统，考虑方程
\[
M \ddot{X} + KX = 0
\]
我们依旧假定解的形式
\[
X= X_0\sin(\omega t + \varphi )
\]
带入
\[
(K- \lambda_i M)X_{i0}=0
\]
其中 \(\lambda_i = \omega^2_i\)
显然这也是一个类似于特征值的求解问题，我们假定解满足
\[
0 \le \lambda_1\le \lambda_2\le \cdots \le\lambda_n
\]
和对应的特征向量（也就是振型）
\[
\varphi_1, \varphi_2，\cdots， \varphi_n
\]
显然，特征值可能出现重根和零根的情况，我们暂时先不考虑
如果出现重根情况，重根对应的特征矢量不能唯一确定（因为其对应的特征矢量任意线性叠加都是其特征矢量），但是我们可以通过正交性确定对应的特征矢量。
我们依旧能给出刚度矩阵、质量矩阵和振型满足的关系式
\[\begin{align}
&\varphi_j^T K\varphi_j =\begin{cases}0 &i\ne j \\ \overline{K}_i &i=j\end{cases}\\
&\varphi_j^T M\varphi_j =\begin{cases}0 &i\ne j \\ \overline{M}_i &i=j\end{cases}\\
&\omega_i = \frac{\overline{ K}_i}{\overline{M}_i}
\end{align}
\]
证明如下：对于第 \(i\) 个和第 \(j\) 个主振型矢量 \(\Phi_{i}\) 和 \(\Phi_{j}\) 有
\[
\begin{align}
&K\Phi_{i}= \lambda_{i}M\Phi_{i} \\
&K\Phi_{j}= \lambda_{j}M\Phi_{j}
\end{align}
\]
对于上述两个式子，分别左乘 \(\Phi_{j}^{T}\) 和 \(\Phi_{i}^{T}\) 可得
\[
\begin{align}
&\Phi_{j}^{T}K\Phi_{i}= \lambda_{i}\Phi_{j}^{T}M\Phi_{i} \\
&\Phi_{i}^{T}K\Phi_{j}= \lambda_{j}\Phi_{i}^{T}M\Phi_{j}
\end{align}
\]
由于 \(K\) 和 \(M\) 都是对称矩阵，所以有
\[
\Phi_{j}^{T}K \Phi_{i} = \Phi_{i}^{T} K \Phi_{j}
\]
\[
\Phi_{j}^{T} M \Phi_{i} = \Phi_{i}^{T}M \Phi_{j}
\]
将其带入之前的表达式中，两式相减可得
\[
(\lambda_{i}-\lambda_{j})\Phi_{i}^{T} M \Phi_{j}=0
\]
所以可以证得正交性。
对于不同固有频率的主振型关于质量矩阵和刚度矩阵正交，这就是多自由度振动系统的主振型的正交性。
有阻尼的多自由度振动的特殊情况与主振型叠加法
由于振型之间相互正交，所以他们构成了 \(N\) 维向量空间的一组正交基
\[
x(t) = \sum_{i=1}^n \varphi_i q_i(t) = T\{q\}
\]
带入到运动方程中有
\[
M(\sum^{n}_{i=1}\varphi_i \ddot{q}_i(t))+C(\sum_{i=1}^n\varphi_i \dot{q}_i(t))+K(\sum_{i=1}^n \varphi_i q_i(t))=F
\]
我们对方程解耦
什么是耦合呢？他表示一个广义坐标下运动对另一个广义坐标下的运动影响，或者说振动元素的相互作用现象，在系统的振动微分方程中，耦合表现为系统的质量矩阵、刚度矩阵或阻尼矩阵的非对角元素不为零。
对于无阻尼振动，我们可以利用正交性来对方程解耦，也就是将存在耦合的振动转换为不存在耦合的情况。解耦后系统的实际运动表现为多个单自由度系统的振动
但是对于阻尼系统，我们不一定能通过正交性来解耦，因为 \(\sum^{n}_{i=1}\varphi_j^T C\varphi_i \dot{q}(t)\) 不一定是一个对角矩阵。一般有两种处理方法：忽略非对角线元素 \(\sum^{n}_{i=1}\varphi_j^T C\varphi_i \dot{q}(t)\approx \overline{C}_jq_j\) (一般适用于弱阻尼系统， \(\xi\) 不大于 0.2 )或采用比例阻尼 \(C = \alpha M+\beta K\)。
于是我们可以将方程解耦成如下形式
\[
\ddot{q}_j+2\xi_j\omega_j\dot{q}_j+\omega_j^2q_j=f_j
\]
其中 \(f_j=F/M_j\)
我们如何确定 \(\xi_j\) 呢？他是一个无量纲常数，一般数值较小。一种办法是直接将其取为一个常数，而对于比例阻尼则有
\[
\xi_j = \frac{1}{2}(\frac{\alpha}{\omega_j}+\beta \omega_j)
\]
我们也可以对初始条件进行解耦，比如说对于位移初始条件：
\[
x(0) = \sum_{i=1}^{n} \varphi_iq_i(0)
\]
左乘 \(\varphi^T_jM\) 有
\[
\varphi_j^TM x(0) = \sum_{i=1}^{n}\varphi_j^TM \varphi_iq_i(0) =\overline{M}_{j}q_{j}(0)
\]
以上方法称为主坐标分析法or主振型叠加法
最后，我们利用杜哈梅积分对上述方程求解
复数解法
考虑受到简谐广义力激励的多自由度系统：
\[
M \ddot{X} + C \dot{X} + KX = F_oe^{i \omega t}
\]
我们可以取
\[
X = X_0 e^{i \omega t}
\]
然后带入到方程求解。
\[
(K-\omega^2M+i\omega C)X_{0}= F_0
\]
解出 \(X = HF_0\)，\(H\) 为复频响应矩阵
\[
H= (K-\omega^2M +i \omega C)^{-1}
\]
所以得到系统响应
\[
x = HF_0e^{i\omega t}
\]
但是这个方法仅适用于小自由度系统（因为要对矩阵求逆，矩阵求逆不是一个很舒服的事情）
复振型解法
Notes
如果要处理带阻尼的多自由度系统就是用这种方法
在引入复振型之前，我们对阻尼矩阵进行一个说明，我们之前引入了两种方法来处理阻尼矩阵：一是只取对角元元素；二是化成比例阻尼处理，但这两个只是极少数的情况，于是在其他情况下，我们至少能保证阻尼矩阵是一个满秩矩阵。（虽然我不知道有什么用哈哈哈）
当存在阻尼作用时，系统一定不是一个保守状态，借用常微分的知识，我们引用一个衰减项
\[
x = X_0e^{\lambda t}
\]
对于自由振动有
\[
[\lambda^2 M+\lambda C+K] e^{\lambda t}=0
\]
要求有解，方程要求
\[
\det[\lambda^2M+\lambda C+K] = 0
\]
这个称为线性阻尼系统的特征方程，他是 \(\lambda\) 的 \(2n\) 次代数方程，它的解 \(\lambda_i\) 可以是实根或复根。如果是实根，他一定是负值，对应系统的衰减自由运动，如果是复根，它一定具有复实部，且复特征根一定是共轭成对出现，对应特定频率与减幅率的一种衰减自由振动。记其对应的特征向量为 \(\phi_{i}\)
我们做一个变换，考虑取
\[
{y}=\left\{\begin{matrix}
\dot{x}\\x
\end{matrix}\right\}
\]
所以方程变成
\[
\widehat{M}\dot{y} +\widehat{K}y = \widehat{F}
\]
其中
\[
\widehat{M} = \begin{bmatrix}
0&M\\
M&C
\end{bmatrix}\quad
\widehat{K}= \begin{bmatrix}
-M&0\\
0&K
\end{bmatrix}\quad
\widehat{F}= \begin{bmatrix}
0\\F
\end{bmatrix}
\]
所以可以设特征振动为
\[
y = \psi e^{\lambda t}
\]
因为我们描述的是同一个运动方程，所以他们有相同的特征根，即
\[
\psi_i=\left\{\begin{matrix}
\lambda_i \phi_i\\\phi_i
\end{matrix}\right\}
\]
我们依旧能证明特征根满足正交性
\[
\begin{align}
&\psi_j^T \widehat{K}\psi_j =\begin{cases}0 &i\ne j \\ \overline{K}_i &i=j\end{cases}\\
&\psi_j^T \widehat{M}\psi_j =\begin{cases}0 &i\ne j \\ \overline{M}_i &i=j\end{cases}\\
&-\lambda_i = \frac{\overline{ K}_i}{\overline{M}_i}
\end{align}
\]
利用正交型和复振型矩阵 \(\Psi = [\psi_1,\dots,\psi_{2n}]\) 设
\[
y = \Psi z
\]
带入方程并左乘 \(\Psi^T\) 解耦得
\[
\Psi^T \widehat{M}\Psi \dot{z} + \Psi^T \widehat{K} \Psi z = \Psi^T \widehat{F}
\]
也就是
\[
\overline{M}_i \dot{z}_i+\overline{K}_iz_i = (\Psi^T \widehat{F})_i
\]
考虑到 \(\overline{K}_{i} = - \lambda_{i} \overline{M}_{i}\) 那么上述方程又可以写作
\[
\dot{z}_i-\lambda_iz_i = \frac{\phi^{T}_{i}F(t)}{\overline{ M } _{i}}
\]
所以对应零初始条件解：
\[
z_{i } = \frac{1}{\overline{ M } _{i}}\int_{0}^{t} \phi^{T}_{i}F(t)e^{\lambda_{i}(t-\tau)} \, {\rm d} \tau
\]
以此回带便可以解出复振型解。
多自由度振动的近似解法
邓克利法
对于保守系统方程，我们可以改写成如下形式
\[
D \ddot{x}+x =0
\]
其中 \(D=FM,F=K^{-1}\) 称为柔度矩阵，我们将 \(x=X_0 e^{i \omega t}\) 带入并改写方程可得
\[
(D-\nu I) X_0 =0
\]
其中 \(\nu= \frac{1}{\omega^2}\) 同样的为了使方程有解我们可以得出
\[
\det [D-\nu I] =0
\]
展开得
\[
\nu^n+\alpha_1\nu^{n-1}+\dots+\alpha_{n-1}+a_n=0
\]
根据特征值的性质，我们知道
\[
\alpha_1 = -(D_{11}+D_{22}+\dots+D_{nn}) = - \text{tr} D
\]
根据多项式的性质，我们可以将特征方程改成如下形式
\[
(\nu-\nu_1)(\nu-\nu_2)\dots(\nu-\nu_n)=0
\]
所以我们有
\[
\sum_{i=1}^{n} \nu_i= \sum^{i=1}_{n}D_{ii}
\]
假定 \(\nu_1\) 为最大值(对应频率最小)，假定其他频率对应的 \(\nu_2,\dots,\nu_n\) 远小于 \(\nu_1\)，因此可以近似忽略，所以我们得到了一个基频近似公式（下限）
\[
\nu_1 \approx \sum_{i=1}^{n}D_{ii}
\]
我们看柔度矩阵 \(F\)，其元素 \(F_{ij}\) 表示 \(i\) 坐标上的单位力引起的质量 \(j\) 的位移。
接下来讨论 \(D\) 。对于 \(D\) 上的对角元素 \(D_{ii} = F_{ii}M_{i}\) 它表示系统内仅有 \(i\) 质量，他的柔度就等于 \(F_{ii}\)，也就是说 系统内仅有第 \(i\) 质量时其固有频率 \(\overset{\sim}{\omega}_i\)满足
\[
D_{ii}=F_{ii}M_{i} = \frac{1}{\overset{\sim}{\omega
}^2_i}
\]
我们就算出了近似基频的取值。
Rayleigh 法
第一 Rayleigh 商：
\[
R_I (x) = \frac{X^T KX}{X^TMX}
\]
第二 Rayleigh 商：
\[
R_{II}(x) = \frac{X^T MX}{X^TMK^{-1}MX}
\]
我们给出 Rayleigh 商满足的关系式
\[
R_{I}(x) \ge R_{II}(x)\ge\omega^2_1
\]
Rayleigh 法给出了最小固有频率的上限，因此结合[[多自由度振动 II ; 多自由度振动的近似计算方法 I#邓利克法|邓利克法]]可以给出最小固有频率的范围。
现在来证明上述不等式，由于一个振型矢量都可以表示为各个特征矢量的线性组合，即：
\[
X=\sum c_i \varphi_i
\]
假定这些特征矢量均已标准化，也就是
\[
\begin{align}
&\varphi^T_i M \varphi_i =1\\
&\varphi^T_iK\varphi_i = \omega_i^2
\end{align}
\]
所以将第一 Rayleigh 商展开得
\[
R_{I}=\frac{(\sum c_i \varphi_i)^TK(\sum c_i \varphi_i)}{(\sum c_i \varphi_i)^TM(\sum c_i \varphi_i)}=\frac{\sum\sum c_ic_j \varphi^T_iK\varphi_j}{\sum\sum c_ic_j \varphi_i^T M \varphi_j}=\frac{\sum c^2_i \omega^2_i}{\sum c_i^2}
\]
我们提取 \(\omega_1\) (最小固有频率)
\[
\frac{\sum c^2_i \omega^2_i}{\sum c_i^2} = \frac{{\sum_{i=2}} c_i (\frac{\omega_i}{\omega_1})^2+c_1}{\sum c^2_i}\omega_1^2\ge \omega_1^2
\]
对于第二 Rayleigh 商，我们先导出特征根关于 \(MK^{-1}M\) 的正交性
\[
\begin{align}
\varphi_j^TMK^{-1}M \varphi_i &= \varphi_j^TM \varphi_i \varphi^{-1}_iK^{-1}(\varphi^T_j)^{-1} \varphi_j^T M \varphi_i\\
& =\begin{cases}
0 &i \ne j\\
\frac{1}{\omega^2_i}&i = j
\end{cases}
\end{align}
\]
简化过程直接给出第二Rayleigh商
\[
\frac{\sum c^2_i}{\sum c^2_i(\frac{1}{\omega_i})^2} = \frac{\sum c^2_i}{\sum_{i=2} c^2_i(\frac{\omega_1}{\omega_i})^2+c^2_1 }\omega^2_1\ge \omega_1^2
\]
下面证明
\[
\frac{\sum c^2_i \omega^2_i}{\sum c_i^2} \ge \frac{\sum c^2_i}{\sum c_i(\frac{1}{\omega_i})^2}
\]
做变换
\[
\sum c_i^2\omega_i^2 \sum c_i(\frac{1}{\omega_i})^2\ge \sum c_i^2\sum c_i^2
\]
由于下标不影响结果，所以展开得
\[
\frac{1}{2}\sum\sum c_i^2 c_j^2 [(\frac{\omega_i}{\omega_j})^2+(\frac{\omega_i}{\omega_j})^2] \ge \sum\sum c_i^2 c_j^2
\]
利用均值不等式即可证明。
上面介绍了两种 Rayleigh 商，前者适用于刚度矩阵已知的情形，后者适用于柔度矩阵已知的情形，一般来说，前者较为简便，后者较为准确。
我们计算出来的 Rayleigh 商虽然不是系统的任一阶固有频率的平方，但必介于系统的最低和最高固有频率的平方 \(\omega_1^2\) 和 \(\omega^2_n\) 之间，即
\[
\omega^2_1 \le R(X) \le \omega^2_n
\]
可以证明，Rayleigh 商在各阶真实振型 \(\varphi_i\) 处取驻点。我们在计算中使用的假设振型越接近系统的真实振型，算出来的固有频率越准确，从物理上理解就是，假设振型相当于对实际系统增加了约束，使系统的刚度提高，因此基频也提高。
Ritz 法
我们将振型表示成有限个独立假设模态的线性和，即
\[
X = \sum_{i=1}^m C_i \psi_i \equiv \Psi \{C\}
\]
其中 \(\psi_i\) 为线性独立的假设模态，\(\Psi\) 为 \(n\times m\) 的矩阵，其中 \(n\) 表示系统内自由度个数，\(m\) 表示假设模态的数量，\(C_i\) 为系数
带到Rayleigh商中
\[
R_{I}= \frac{\{C\}^T\Psi^T K \Psi \{C\}}{\{C\}^T \Psi^T M \Psi \{C\}}=\frac{\{C\}^T \overline{K}\{C\}}{\{C\} \overline{M} \{C\}} = \lambda
\]
我们要使我们的 Rayleigh 商最小，即取驻点 \(\frac{\partial R_{I}}{\partial {C}_i}=0\)，我们可以得出关系式
\[
[\overline{K}-\lambda \overline{M}] \{C\}=0
\]
可求出对应的特征值 \(\lambda_1\dots \lambda_m\) 和特征向量 \(C_1\dots C_m\)
我们同样可以用 Rayleigh 第二商求解。
事实上我们所得出的解坐落在一个线性空间 \(\text{span}\{ \varphi_1,\varphi_2 \dots \varphi_m\}\) 中，而这个空间是原本 \(n\) 维空间的子空间，Ritz 法实际上起着坐标缩并的作用。用 Ritz 法计算处的固有频率前半频率精度较高，因此要计算前 \(s\) 个固有频率，通常取 \(r=2s\) 个假设振型。
矩阵迭代法
将系统的主振型方程可表示为柔度形式，即
\[
DX_i = \nu_i X_i
\]
所以我们可采用如下迭代算法更新
\[
\begin{align}
1.&\text{任取振型} x_0\\
2.&第k次迭代:\\
& \beta _kx_k = Dx_{k-1}\\
\end{align}
\]
\(\beta_k\) 表示向量 \(Dx_{k-1}\) 迭代后所产生向量中的最大分量，其中 \(x_k\) 满足最大分量为1
所以我们可以的得到近似关系
\[
\begin{align}
&\beta \rightarrow \nu_1 = \frac{1}{\omega_1^2}\\
&x_n \rightarrow \varphi_1
\end{align}
\]
也就是计算方法中的改进幂法。
如果我们用刚度矩阵方式进行迭代即方程
\[
M^{-1}K X_i = \lambda_iX_i
\]
所求得的固有频率为最大固有频率。
如果我们想求解 \(\nu_2,\varphi_2\) 则需要清除假设振型中 \(\varphi_1\) 的影响。 对于任意振型
\[
X = a_1 \varphi_1+\sum_{i=2} a_i \varphi_i
\]
利用正交性，左乘 \(\varphi_1^T M\) 得
\[
\varphi^T_1 M X = \varphi^T_1 M(a_1\varphi_1+\sum_{i=2} a_i \varphi_i) = a_1 \varphi_1^T M \varphi_1
\]
考虑消除 \(\varphi_1\) 得振型
\[
X' = X- a_1 \varphi_1
\]
由于 \(a_1\) 是一个数，所以 \(\varphi_1\) 的位置并不会影响，我们将正交性结果代入得
\[
X' = X- \varphi_1 \frac{\varphi_1^TMX}{\varphi^T_1M \varphi_1 } =(I- \frac{\varphi_1 \varphi^T_1M }{\varphi_1^T M \varphi_1}) X
\]
所以我们可以得到清形矩阵
\[
H= I - \frac{\varphi_1 \varphi^T_1M }{\varphi_1^T M \varphi_1}
\]
推广到 \(r\) 个自由度为
\[
H_r= I - \sum_{i=1}^r\frac{\varphi_i \varphi^T_iM }{\varphi_i^T M \varphi_i}
\]
子空间迭代法
对于一个具有 \(n\) 个自由度的多自由振动系统，假定我们需要求取 \(l\) 个振型，我们一般会取 \(m\) 个振型，进行如下迭代：
先选取任意基向量
\[
\Psi_{n\times m} = [\varphi_1,\varphi_2,\dots,\varphi_m]
\]
左乘矩阵 \(D\)
\[
\Psi_1 = D \Psi_0
\]
此时我们进行正交化处理：
\[
\begin{align}
&\overline{K} = \Psi_1^TK \Psi_1\\
&\overline{M}= \Psi_1^T M \Psi_1
\end{align}
\]
按照之前的方法，我们可以得到此时对应的 \(m\) 个特征值和特征值对应的特征向量
\[
\begin{align}
&\lambda_1,\lambda_2\dots,\lambda_m\\
&\overline{C}_1,\overline{C}_2\dots,\overline{C}_m
\end{align}
\]
我们按照如下方法构造新的基向量
\[
\Psi^*_{1} = [\Psi_1]\begin{bmatrix}
\overline{C}_1&\overline{C}_2&\dots&\overline{C}_m
\end{bmatrix}
\]
因此我们将上述过程简化成如下
\[
\begin{align}
1.&任意选取初始基向量 \Psi^*_0\\
2.&进行如下迭代:\\
&a. 计算 \Psi_{i} = D \Psi^*_{i-1}\\
&b.计算系数矩阵 C = \begin{bmatrix}
\overline{C}_1&\overline{C}_2&\dots&\overline{C}_m
\end{bmatrix}\\
&c.计算\Psi^*_i=\Psi_iC
\end{align}
\]
当前 \(l\) 个特征值(也就是我们所需要的)满足如下收敛条件
\[
\left|{\lambda^{(k)}_i-\lambda_i^{(k-1)}\over \lambda^{(k)}_i}\right| \le \varepsilon
\]
时，我们判断已成功求出我们所想要的振型。一般而言 \(\varepsilon\) 取 \(10^{-3}\) 或 \(10^{-4}\) 即可满足条件，\(\varepsilon\) 越小，所需要的计算量级越大。
这种方法称作是子空间迭代法。
也就是说，我们在迭代过程中求得 \(\Psi_1\) 后，并不直接对它进行迭代，而是在迭代前先对 \(\Psi_1\) 进行处理。首先先进行一个正交化，这样可以使它的各列经迭代后分别趋于各个不同阶的主振型，而不是都趋于一阶主振型，其次，还要对它进行标准化或基准化，使得在数字运算中能保持适当的大小。
我们知道，子空间迭代法效率与 \(m\) 的选取有关，一般而言 \(m\) 越大，振型计算结果越精确，但其所需要的计算量也越大；而 \(m\) 越小，计算量越小，但收敛速度越慢。我们一般使 \(m\) 和 \(l\) 满足如下关系
\[
m =\min\{2l,l+8\}
\]
在计算软件，我们采用如下办法确定我们的 \(\Psi_0\)。我们尽可能确保全部信息，所以取 \(\varphi_1 = \begin{bmatrix}1&1&\dots&1\end{bmatrix}^T\)，而对剩下的分量，我们取质量矩阵和刚度矩阵的对角元素并做除法得 \(K_{ii}\over M_{ii}\) 按照其大小从小到大排序(因为基频越低越是我们想要的元素)，按照大小及其对应的位置确定分量：比如说最小的 \(K_{ii}\over M_{ii}\) 对应着第 5 个元素（也就是 \(i=5\) ），\(\varphi = \begin{bmatrix}0&0&0&0&1&0&\dots&0\end{bmatrix}\) (真的是很神奇的确定方法)
对于非正定或者是病态的刚度矩阵 \(K\) ，我们一般做一个移位处理，即
\[
[K-\omega^2M]X = [(K+\alpha M)-(\omega^2+\alpha)M]X = 0
\]
这样我们就得到一个等效的刚度矩阵 \(\overline{K} = K+\alpha M\) 这个矩阵显然是对称且正定的。我们需要合理的选取 \(\alpha\) ，如果 \(\alpha\) 过大就会减小刚度矩阵的影响，也就是 \(M\) 吃掉了 \(K\) ；如果 \(\alpha\) 过小，显然这有可能依旧是一个病态矩阵，我们一般根据对角元素来确定合适的 \(\alpha\)。此外，合适的 \(\alpha\) 有利于我们提高收敛速度。