如何用Python验证Cayley-Hamilton定理？NumPy实战教程

最新推荐文章于 2026-07-01 17:22:11 发布

原创

最新推荐文章于 2026-07-01 17:22:11 发布 · 719 阅读

用Python亲手“触摸”Cayley-Hamilton定理：从抽象公式到NumPy代码的奇妙旅程

你是否曾觉得线性代数里的那些定理，比如Cayley-Hamilton定理，听起来高深莫测，仿佛悬浮在纯数学的云端，与你的代码世界隔着一条鸿沟？我记得自己第一次看到“一个矩阵代入其自身的特征多项式结果为零矩阵”这个表述时，也是一头雾水。直到有一天，我决定不再仅仅阅读证明，而是打开Python，用NumPy亲手去“计算”这个定理。那一刻，抽象的数学符号瞬间变成了屏幕上可验证的、实实在在的数字。这不仅让我彻底理解了定理的内涵，更让我体会到，将数学理论与编程实践相结合，是掌握复杂概念最有力的方式。这篇文章，就是为你——一位对数学好奇、并习惯用代码探索世界的Python开发者——准备的一次动手实验。我们将不满足于理论推导，而是共同编写代码，见证一个矩阵如何“吞噬”掉自己的特征多项式，最终归于零矩阵的奇妙景象。

1. 出发之前：理解我们即将验证的核心

在开始敲代码之前，我们有必要用更“程序员友好”的语言，重新梳理一下Cayley-Hamilton定理到底在说什么。这能确保我们的代码验证不是盲目的计算，而是有目的的探索。

简单来说，对于任何一个 n×n 的方阵 A，我们都可以为它定义一个独一无二的“身份证”——特征多项式。这个多项式通常写作：

χ_A(λ) = det(λI - A) = λ^n + c_{n-1}λ^{n-1} + ... + c_1λ + c_0

这里，det 表示行列式，I 是单位矩阵，λ 是一个变量（通常用 s 或 λ 表示），而 c_i 则是计算出来的系数。

Cayley-Hamilton定理则提出了一个非常反直觉的断言：如果你把这个多项式里的变量 λ，直接替换成矩阵 A 本身，然后按照矩阵运算的规则去计算这个多项式，那么最终得到的结果，将是一个零矩阵。

用公式表示就是： χ_A(A) = A^n + c_{n-1}A^{n-1} + ... + c_1A + c_0I = 0 (这里的 0 是零矩阵)

注意：这里的 A^n 是矩阵的 n 次幂，c_0I 表示常数项 c_0 乘以单位矩阵。这是矩阵多项式计算的关键，标量运算和矩阵运算必须区分清楚。

这个定理的奇妙之处在于，它揭示了矩阵自身蕴含的一种强烈的“自洽性”。一个矩阵居然能满足一个以自己为系数的多项式方程。从应用角度看，这个定理为简化矩阵的高次幂计算、求解矩阵函数（如 e^A, sin(A)）以及系统控制理论中的分析提供了重要的理论工具。

为了在Python中验证它，我们的任务清晰地分为三步：

计算特征多项式：给定矩阵 A，找出它的特征多项式系数 c_i。
构造矩阵多项式：用矩阵 A 和计算出的系数，按照定理公式计算 χ_A(A)。
验证零矩阵：检查计算结果是否足够接近零矩阵（由于浮点数计算，我们检查其元素是否都接近0）。

2. 搭建实验环境与核心工具函数

工欲善其事，必先利其器。我们首先确保环境就绪，并构建几个核心函数，它们将是我们验证过程的基石。

2.1 环境准备与基础导入

我们将完全依赖 NumPy 和 SciPy 这两个科学计算的核心库。SciPy 的 linalg 模块提供了更丰富的线性代数函数。

import numpy as np
from scipy import linalg
import warnings
warnings.filterwarnings('ignore')  # 暂时忽略一些不影响计算的警告

2.2 核心函数一：获取特征多项式系数

计算特征多项式系数有多种方法。我们将实现两种最直接的方法，并对比其异同，这本身也是一个很好的学习过程。

方法A：通过特征值构建多项式 理论上，特征多项式的根就是矩阵的特征值。如果矩阵 A 的特征值为 λ1, λ2, ..., λn，那么其特征多项式为 (λ - λ1)(λ - λ2)...(λ - λn)。展开后即可得到系数。但需要注意的是，通过根反推系数在数值计算上可能对舍入误差更敏感。

def char_poly_from_eigvals(A):
    """
    通过计算特征值来推导特征多项式系数。
    注意：对于具有重复特征值或病态矩阵，此方法数值稳定性可能较差。
    """
    eigvals = np.linalg.eigvals(A)
    # np.poly 函数根据给定的根返回多项式系数，最高次项系数为1
    coefficients = np.poly(eigvals)
    return coefficients

方法B：使用Faddeev-Leverrier算法 这是一种经典的、通过矩阵的迹来递推计算特征多项式系数的数值算法。它不直接计算特征值，通常具有更好的数值性质。

def char_poly_faddeev_leverrier(A):
    """
    使用 Faddeev-Leverrier 算法计算特征多项式系数。
    返回系数列表 [c_n, c_{n-1}, ..., c_0]，其中 c_n = 1。
    """
    n = A.shape[0]
    coeffs = np.ones(n + 1)  # 初始化系数数组，最高次项系数为1
    M = A.copy()
    for k in range(1, n + 1):
        # 计算当前矩阵的迹，并除以k得到系数
        ck = -M.trace() / k
        coeffs[n - k] = ck
        if k < n:
            M = A @ (M + ck * np.eye(n))
    # 最后一项是常数项
    coeffs[-1] = -M.trace() / n if n > 0 else 1
    # 调整符号，使得特征多项式为 det(λI - A) = λ^n + c_{n-1}λ^{n-1}+...+c_0
    # Faddeev-Leverrier 通常给出 det(A - λI) 的系数，所以需要调整
    coeffs = (-1)**np.arange(n+1) * coeffs
    return coeffs

提示：np.poly 函数虽然方便，但它是通过构造伴随矩阵来求特征值，对于稍大的矩阵或条件数大的矩阵，误差可能累积。Faddeev-Leverrier 算法是一种直接的方法，值得掌握。

为了更直观地对比这两种方法，我们设计一个简单的测试：

def test_coeff_methods():
    # 用一个简单的 3x3 矩阵测试
    A_test = np.array([[2, -1, 0],
                       [-1, 2, -1],
                       [0, -1, 2]])
    coeffs_eig = char_poly_from_eigvals(A_test)
    coeffs_fl = char_poly_faddeev_leverrier(A

标签