Eigen — 特征值与模态分析¶

概述¶

结构动力学的特征值计算和模态分析。TinyMath 依赖链中的第 4 层。

头文件： #include "tiny_eigen.h" / #include "tiny_modal.h"

架构¶

该模块分为两个互补部分：

子模块	头文件	用途
特征值	`tiny_eigen.h`	2×2 闭式解、Schur 提取、Hessenberg 约化、Francis 双移位 QR
模态	`tiny_modal.h`	MIMO 模态分析：极点 → 频率/阻尼、MAC、模态滤波/排序

处理流程¶

A (n×n) → Hessenberg → QR 迭代 → Schur 形式 → 特征值 (tiny_cfloat_t[])
                                        ↓
                              模态分析：
                              极点 → ω_n, ζ → MAC → 滤波/排序

特征值计算¶

2×2 闭式解通用矩阵 (\(n \times n\))Schur 提取

tiny_error_t tiny_eig_2x2_f32(float a00, float a01, float a10, float a11,
                               tiny_cfloat_t *lam1, tiny_cfloat_t *lam2);

\(2 \times 2\) 矩阵的直接解析解。无需迭代。

参数	说明
`a00, a01, a10, a11`	按行优先顺序的矩阵元素
`lam1, lam2`	输出特征值（复数）

公式： \(\lambda = \dfrac{\text{tr}(A) \pm \sqrt{\text{tr}(A)^2 - 4\det(A)}}{2}\)

\(2 \times 2\) 子块最快路径。Francis QR 迭代内部用于 \(2 \times 2\) Wilkinson 移位。

size_t       tiny_eig_general_f32_workspace_size(int n);
tiny_error_t tiny_eig_general_f32(const float *A, int n,
                                   tiny_cfloat_t *eigs,
                                   void *ws, size_t ws_size);

处理流程： \(A \xrightarrow{\text{Hessenberg}} H \xrightarrow{\text{Francis QR}} \text{Schur} \xrightarrow{\text{提取}} \Lambda\)

参数	说明
`A`	输入 \(n \times n\) 矩阵（不修改）
`n`	矩阵维度
`eigs`	输出 \(n\) 个复数特征值数组
`ws`	工作空间缓冲区（\(\geq\) `tiny_eig_general_f32_workspace_size(n)` 字节）
`ws_size`	工作空间大小（安全校验）

tiny_error_t tiny_schur_extract_eigs_f32(const float *T, int n, int step,
                                          tiny_cfloat_t *eigs);

从 拟三角（实 Schur）矩阵 \(T\) 中提取特征值。

参数	说明
`T`	\(n \times n\) 拟三角矩阵，行优先
`step`	行步长（列数 + padding）
`eigs`	输出特征值

实特征值在对角线上；复共轭对出现在 \(2 \times 2\) 块中。

模态分析¶

极点 → 频率和阻尼MAC（模态保证准则）模态滤波

tiny_error_t tiny_modal_poles_to_modal_f32(const tiny_cfloat_t *poles,
                                            float *freq, float *damping,
                                            int n, float fs);

将离散时间极点转换为固有频率和阻尼比：

参数	说明
`poles`	\(n\) 个复极点（来自 SSI/ERA）
`freq`	输出：固有频率（Hz）
`damping`	输出：阻尼比（%）
`n`	极点数量
`fs`	采样频率（Hz）

公式： \(\(p = \ln(\lambda) \cdot f_s, \quad \omega_n = |p|, \quad \zeta = -\frac{\Re(p)}{|p|}\)\)

tiny_error_t tiny_modal_mac_matrix_f32(const float *Phi1, int n1,
                                        const float *Phi2, int n2,
                                        int n_dof, float *MAC);
float tiny_modal_mac_pair_f32(const float *phi_a, const float *phi_b, int n_dof);

函数	说明
`mac_matrix`	两组振型之间的完整 \(n_1 \times n_2\) MAC 矩阵
`mac_pair`	两个振型向量之间的单一 MAC 值

公式： \(\text{MAC}(\phi_a, \phi_b) = \dfrac{|\phi_a^T \phi_b|^2}{(\phi_a^T \phi_a)(\phi_b^T \phi_b)}\)

MAC = 1 表示完全相关模态；MAC = 0 表示正交模态。

int tiny_modal_filter_stable_poles(const tiny_cfloat_t *poles,
                                    int n, float max_damping,
                                    float min_freq, float max_freq,
                                    float fs);
int tiny_modal_sort_by_freq(const tiny_cfloat_t *poles,
                             const float *freq, const float *damping,
                             int n, int *order);

函数	说明
`filter_stable_poles`	返回满足 \(\zeta < \text{max\_damping}\) 且 \(\omega \in [\text{min\_freq}, \text{max\_freq}]\) 的极点数
`sort_by_freq`	返回按固有频率升序排列的置换数组

典型 SSI/ERA 工作流

// 1. 从系统矩阵计算特征值
tiny_cfloat_t poles[50];
tiny_eig_general_f32(A, 50, poles, ws, ws_size);

// 2. 转换为频率和阻尼
float freq[50], damping[50];
tiny_modal_poles_to_modal_f32(poles, freq, damping, 50, fs);

// 3. 滤波稳定极点
int n_stable = tiny_modal_filter_stable_poles(poles, 50, 5.0f, 0.1f, 50.0f, fs);

// 4. 计算两次识别结果的 MAC
float mac_matrix[25 * 25];
tiny_modal_mac_matrix_f32(Phi_run1, 25, Phi_run2, 25, n_dof, mac_matrix);

算法细节¶

Francis 双移位 QR¶

\(n \times n\) 特征值计算的核心算法：

Hessenberg 约化： 通过 Householder 反射将 \(A \to H\)（\(O(n^3)\)）
Francis 双移位 QR： 在 \(H\) 上隐式双移位迭代（每次 \(O(n^3)\)）
收缩： 当 \(H_{n,n-1} \approx 0\) 时，分割为 \((n-1) \times (n-1)\) 和 \(2 \times 2\)
Schur 提取： 从最终的拟三角矩阵中提取特征值

双移位保留了实矩阵的实运算——复特征值以 \(2 \times 2\) 块的形式出现在对角线上。

收敛性¶

条件良好的矩阵通常 2-3 次 QR 迭代收敛。特征值聚集的矩阵可能需要更多迭代。算法使用 Wilkinson 移位 加速收敛。