Eigen — 实现¶

文件结构¶

eigen/
├── tiny_eigen.h        (15 行 — 特征值 API)
├── tiny_eigen.c        (210 行 — 特征值实现)
├── tiny_modal.h        (24 行 — 模态分析 API)
└── tiny_modal.c        (157 行 — 模态分析实现)

依赖：tiny_cfloat.h、tiny_decomp.h、<math.h>。

特征值处理流程¶

2×2 闭式解¶

最简单的路径——无迭代，无需工作空间：

float tr  = a00 + a11;
float det = (a00 - a11)*(a00 - a11) + 4*a01*a10;

if (det >= 0) {              // 实特征值
    lam = (tr ± sqrt(det)) / 2;
} else {                     // 复共轭对
    re = tr / 2;
    im = sqrt(-det) / 2;
}

作为 Francis QR 迭代中收缩的基础情况使用。

Schur 提取¶

Francis QR 迭代后，\(A\) 成为 实 Schur 形式（拟三角）：

实 Schur (n=4)：
┌───┬───┬───┬───┐
│ x │ x │ x │ x │
│ 0 │ x │ x │ x │  ← 次对角线元素为 0（或接近 0）
│ 0 │ 0 │ x │ x │
│ 0 │ 0 │ x │ x │  ← 2×2 块 → 复共轭对
└───┴───┴───┴───┘

函数遍历对角线：

• 若 T[k+1][k] ≈ 0 → 在 T[k][k] 处为实特征值
• 若 T[k+1][k] ≠ 0 → \(2 \times 2\) 块 → 通过 tiny_eig_2x2_f32 提取复共轭对

通用特征分解¶

tiny_eig_general_f32(A, n, eigs, ws, ws_size)

算法流程：

1. 分配工作空间（n×n 矩阵 + n 向量 + n² 暂存）
2. 复制 A → H（行优先，连续）
3. Hessenberg 约化 H（Householder 反射）
4. Francis 双移位 QR 迭代 H（隐式移位，Wilkinson 移位）
5. 从 H（拟三角）中提取特征值
6. 返回特征值到 eigs[]

工作空间布局：

模态分析¶

极点转换¶

// 从离散特征值 λ 到连续时间极点：
p = ln(λ) * fs

// 固有频率和阻尼：
ω_n = |p|              // rad/s
f  = ω_n / (2π)        // Hz
ζ  = -Re(p) / |p|      // 阻尼比（0 = 无阻尼，1 = 临界阻尼）

MAC 计算¶

MAC(φ_a, φ_b) = |φ_a^T φ_b|² / ((φ_a^T φ_a)(φ_b^T φ_b))

• 配对： 两个向量的直接点积，\(O(n_{\text{dof}})\)
• 矩阵： \(n_1 \times n_2\) 两两组合，\(O(n_1 \cdot n_2 \cdot n_{\text{dof}})\)

稳定极点滤波¶

极点满足以下条件时被视为"稳定"：

函数统计满足条件的极点数——调用者据此提取对应的子集。

数值注意事项¶

• Francis QR 使用 Wilkinson 移位 加速收敛——移位从右下角 \(2 \times 2\) 子矩阵计算
• 收缩阈值： \(10^{-14} \cdot \max(1, |T_{kk}| + |T_{k+1,k+1}|)\)
• 仅实矩阵： 算法保留实运算——提取特征值前无复数开销
• Float32 精度限典型 SSI 矩阵（条件数 \(\sim 10^6\)）的特征值精度约为 \(\sim 10^{-4}\)