用Python和NumPy手把手教你计算广义特征向量（附完整代码与可视化）

在工程计算和数据分析中，我们经常会遇到无法对角化的"缺陷矩阵"。这时，广义特征向量就成为了解决问题的关键钥匙。本文将带你用Python的NumPy和SciPy库，从基础概念到代码实现，一步步掌握这个强大的数学工具。

1. 理解广义特征向量的核心概念

1.1 从标准特征向量到广义特征向量

标准特征向量的定义我们都熟悉：对于矩阵A和标量λ，若存在非零向量v使得Av = λv，则v称为A的特征向量，λ为对应的特征值。这个定义可以等价表示为：

import numpy as np

A = np.array([[2, 1, -1], 
              [1, 2, -1], 
              [-1, -1, 2]])
λ = 1
v = np.array([1, 0, 1])  # 示例特征向量

# 验证特征向量
print(np.allclose(A @ v, λ * v))  # 输出True

但当矩阵A不能对角化时，标准特征向量就不足以描述完整的变换行为。这时我们需要 广义特征向量 ——满足(A-λI)ʲv=0的最小j对应的向量v。

1.2 广义特征空间的分层结构

广义特征空间呈现出有趣的层级关系：

层级j	空间Eⱼ(λ)	描述
1	Ker(A-λI)	标准特征空间
2	Ker((A-λI)²)	包含一阶广义特征向量
...	...	...
k	Ker((A-λI)ᵏ)	广义特征空间稳定

这种层级结构在解决微分方程和控制理论问题时特别有用。

2. 计算广义特征向量的Python实现

2.1 准备工作：安装必要库

确保你的Python环境已安装以下库：

pip install numpy scipy matplotlib

2.2 分步计算过程

让我们通过一个具体例子演示计算过程。考虑矩阵：

A = np.array([[3, 1, 1],
              [-2, 0, -2],
              [1, 1, 3]])
λ = 2  # 已知特征值

步骤1：计算标准特征空间

I = np.eye(3)
B = A - λ * I
# 计算零空间
from scipy.linalg import null_space
standard_eigenvectors = null_space(B)
print("标准特征向量：\n", standard_eigenvectors)

步骤2：寻找广义特征向量

当标准特征向量不足时，我们需要计算更高阶的零空间：

B_power = np.linalg.matrix_power(B, 2)
generalized_eigenvectors = null_space(B_power)
print("广义特征向量：\n", generalized_eigenvectors)

2.3 自动化计算函数

下面是一个完整的函数，可以计算任意矩阵的广义特征向量：

def compute_generalized_eigenvectors(A, λ, max_iter=10):
    I = np.eye(A.shape[0])
    B = A - λ * I
    prev_nullity = 0
    eigenvectors = []
    
    for j in range(1, max_iter+1):
        B_j = np.linalg.matrix_power(B, j)
        current_eigenvectors = null_space(B_j)
        current_nullity = current_eigenvectors.shape[1]
        
        if current_nullity > prev_nullity:
            eigenvectors.append((j, current_eigenvectors))
            prev_nullity = current_nullity
        else:
            break
            
    return eigenvectors

3. 可视化特征空间变换

3.1 二维投影可视化

对于三维矩阵，我们可以通过投影观察特征空间的变化：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def plot_eigen_spaces(A, λ):
    eigenvectors = compute_generalized_eigenvectors(A, λ)
    
    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111, projection='3d')
    
    # 绘制单位球面
    u = np.linspace(0, 2 * np.pi, 100)
    v = np.linspace(0, np.pi, 100)
    x = np.outer(np.cos(u), np.sin(v))
    y = np.outer(np.sin(u), np.sin(v))
    z = np.outer(np.ones(np.size(u)), np.cos(v))
    ax.plot_surface(x, y, z, color='b', alpha=0.1)
    
    # 绘制特征向量
    colors = ['r', 'g', 'y']
    for i, (j, vecs) in enumerate(eigenvectors):
        for k in range(vecs.shape[1]):
            v = vecs[:, k]
            ax.quiver(0, 0, 0, v[0], v[1], v[2], 
                     color=colors[i], 
                     label=f'j={j}' if k==0 else "")
    
    ax.set_title(f'特征值λ={λ}的特征空间层级')
    ax.legend()
    plt.show()

plot_eigen_spaces(A, 2)

3.2 特征空间维度变化图

我们可以绘制随着j增大，零空间维度的变化情况：

def plot_nullity_progression(A, λ):
    nullities = []
    max_j = 5
    I = np.eye(A.shape[0])
    B = A - λ * I
    
    for j in range(1, max_j+1):
        B_j = np.linalg.matrix_power(B, j)
        nullity = B_j.shape[0] - np.linalg.matrix_rank(B_j)
        nullities.append(nullity)
    
    plt.figure(figsize=(8, 5))
    plt.plot(range(1, max_j+1), nullities, 'bo-')
    plt.xlabel('j (幂次)')
    plt.ylabel('零空间维度')
    plt.title(f'λ={λ}时(A-λI)^j的零空间维度变化')
    plt.grid(True)
    plt.show()

plot_nullity_progression(A, 2)

4. 工程应用实例：振动系统分析

4.1 缺陷刚度矩阵的处理

在结构动力学中，我们经常会遇到形如Mx'' + Kx = 0的方程。当刚度矩阵K有缺陷时，广义特征向量就派上用场了。

考虑一个简单的质量-弹簧系统：

M = np.diag([2, 2, 2])  # 质量矩阵
K = np.array([[4, -2, 0], 
              [-2, 4, -2], 
              [0, -2, 4]])  # 刚度矩阵

# 求解广义特征问题
from scipy.linalg import eig
λs, vs = eig(K, M)

4.2 构建完整的解空间

当几何重数小于代数重数时，我们需要补充广义特征向量来构建完整的解空间：

def build_complete_solution(A, eigenvalues):
    complete_basis = []
    for λ in eigenvalues:
        # 计算标准特征向量
        B = A - λ * np.eye(A.shape[0])
        std_vecs = null_space(B)
        complete_basis.extend(std_vecs.T)
        
        # 检查是否需要广义特征向量
        if std_vecs.shape[1] < np.sum(np.isclose(eigenvalues, λ)):
            # 计算广义特征向量
            gen_vecs = null_space(np.linalg.matrix_power(B, 2))
            additional_vecs = [v for v in gen_vecs.T 
                             if not any(np.allclose(v, sv) for sv in std_vecs.T)]
            complete_basis.extend(additional_vecs)
    
    return np.vstack(complete_basis).T

4.3 系统响应模拟

利用广义特征向量，我们可以更准确地模拟系统的动态响应：

def simulate_system_response(M, K, initial_conditions, t):
    # 求解广义特征问题
    λs, vs = eig(K, M)
    Φ = build_complete_solution(np.linalg.inv(M) @ K, λs)
    
    # 初始条件变换到模态坐标
    q0 = np.linalg.pinv(Φ) @ initial_conditions[0]
    q0_dot = np.linalg.pinv(Φ) @ initial_conditions[1]
    
    # 计算各模态响应
    response = np.zeros((len(t), M.shape[0]))
    for i in range(Φ.shape[1]):
        ω = np.sqrt(np.abs(λs[i]))
        mode_response = (q0[i] * np.cos(ω * t) + 
                        q0_dot[i]/ω * np.sin(ω * t))
        response += np.outer(mode_response, Φ[:, i])
    
    return response

5. 性能优化与数值稳定性

5.1 处理数值误差的技巧

在实际计算中，我们需要注意：

• 特征值多重性的判断应使用容差而非精确相等
• 零空间计算时设置适当的阈值
• 幂次计算时的数值稳定性

改进后的计算函数：

def robust_null_space(matrix, tol=1e-10):
    u, s, vh = np.linalg.svd(matrix)
    null_mask = (s <= tol)
    return vh[null_mask].T

def stable_matrix_power(B, power):
    # 使用SVD分解提高计算稳定性
    u, s, vh = np.linalg.svd(B)
    return u @ np.diag(s**power) @ vh

5.2 大型稀疏矩阵的处理

对于大型系统，我们可以利用稀疏矩阵特性：

from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigs, svds

def sparse_generalized_eigen(A_sparse, λ, k=3):
    I_sparse = csr_matrix(np.eye(A_sparse.shape[0]))
    B = A_sparse - λ * I_sparse
    _, s, vh = svds(B, k=k, which='SM')
    mask = s < 1e-10
    return vh[mask].T

6. 常见问题与调试技巧

6.1 诊断计算问题

当计算结果不符合预期时，可以检查：

1. 特征多项式是否正确：

   char_poly = np.poly(A)
   print("特征多项式系数：", char_poly)
   

2. 几何重数与代数重数是否匹配：

   def check_multiplicity(A, λ):
       algebraic = np.sum(np.isclose(np.linalg.eigvals(A), λ))
       geometric = A.shape[0] - np.linalg.rank(A - λ*np.eye(A.shape[0]))
       return algebraic, geometric
   

6.2 处理病态矩阵

对于条件数大的矩阵：

def regularized_eigen(A, λ, ε=1e-6):
    I = np.eye(A.shape[0])
    B = A - λ * I
    # 添加小量正则化
    B_reg = B + ε * I
    return null_space(B_reg)

7. 扩展应用：控制系统分析

在控制理论中，广义特征向量用于分析系统可控性和构造状态观测器。例如：

def controllability_matrix(A, B):
    n = A.shape[0]
    C = B
    for i in range(1, n):
        C = np.hstack((C, np.linalg.matrix_power(A, i) @ B))
    return C

def is_controllable(A, B):
    C = controllability_matrix(A, B)
    return np.linalg.matrix_rank(C) == A.shape[0]

8. 进阶主题：广义特征向量的其他应用

8.1 矩阵函数计算

利用广义特征向量可以计算矩阵函数，如矩阵指数：

def matrix_exp_using_generalized_eigen(A, t):
    λs, vs = np.linalg.eig(A)
    Φ = build_complete_solution(A, λs)
    J = np.linalg.inv(Φ) @ A @ Φ  # Jordan形式
    return Φ @ np.diag(np.exp(λs * t)) @ np.linalg.inv(Φ)

8.2 微分方程求解

对于常系数线性微分方程组：

def solve_ode_system(A, x0, t):
    λs, vs = np.linalg.eig(A)
    Φ = build_complete_solution(A, λs)
    c = np.linalg.solve(Φ, x0)
    solution = np.zeros((len(t), len(x0)), dtype=complex)
    for i, λ in enumerate(λs):
        solution += np.outer(c[i] * np.exp(λ * t), Φ[:, i])
    return solution.real