# 教材位置：第8章，第8.2节；图 8-1。
# 单位、变量与形状：A 为 (2,2)，特征向量按列组成 (2,2)；坐标单位统一，特征值是轴向缩放倍
# 数。
# 案例 25：用对称正定矩阵理解特征值、特征向量与几何变换。
# 学习目标：把 Av=λv 的代数关系与“沿特征方向只发生伸缩”的图像联系起来，
# 并检验实对称矩阵的特征向量可以构成一组标准正交基。
# 核心关系：A=Q diag(λ) Qᵀ，tr(A)=∑λ，det(A)=∏λ；本例两个特征值均为正。
# 数据与形状：A 为 2×2 矩阵；vectors 的每一列是一个特征向量；
# 单位圆用一个 2×200 数组表示，每一列都是平面上的一个点。
# 输出：打印特征值、迹和行列式；断言核对恒等式；生成 e25_eigen.png。
# 运行方式：安装本目录 requirements.txt 中的依赖，在本目录执行
# python e25_eigen.py。图片保存到运行时的当前工作目录，不弹出交互窗口。

# NumPy 负责矩阵计算；约定 @ 表示矩阵乘法，* 表示逐元素运算或广播。
import numpy as np
# Matplotlib 负责把单位圆及其线性变换后的椭圆画在同一坐标系内。
import matplotlib.pyplot as plt
# 在创建第一幅图之前配置可跨系统选择的中文字体。
from plot_style import configure_chinese_plot

configure_chinese_plot()
# A=Aᵀ，因此适用实对称矩阵专用的 eigh；浮点数便于后续数值线性代数计算。
A = np.array([[3.0, 1.0], [1.0, 2.0]])
# eigh 按特征值从小到大返回 values，形状为 (2,)；vectors 形状为 (2,2)。
# 第 j 列 vectors[:,j] 对应 values[j]，不是第 j 行；各列长度为 1 且两两正交。
# 特征向量的正负号不唯一，换成相反方向仍是同一个特征方向。
values, vectors = np.linalg.eigh(A)
# 手算 det(A-λI)=λ²-5λ+5=0，得到 (5±√5)/2；按升序排好以与 eigh 对齐。
# 两根均为正，所以 A 在任意非零方向的二次型 xᵀAx 都为正。
manual = np.array([(5 - np.sqrt(5)) / 2, (5 + np.sqrt(5)) / 2])
# 只在显示时保留 6 位小数；后续检验继续使用未舍入的原始数值。
print("特征值 =", np.round(values, 6))
# 迹为主对角线元素之和 5，行列式为 3×2-1×1=5，可独立核对特征值结果。
print("迹、行列式 =", np.trace(A), round(np.linalg.det(A), 6))
# 浮点计算可能有舍入误差，使用 allclose 的容差比较数值，不要求逐位相同。
assert np.allclose(values, manual)
# 同时检验两条 Av_j=λ_jv_j：右侧 (2,2)*(2,) 把第 j 列乘以 λ_j。
# 这条广播规则十分关键；如果误把特征向量当作行，关系一般不成立。
assert np.allclose(A @ vectors, vectors * values)
# QᵀQ=I 检验列向量标准正交；因此在这个实对称案例中 Q⁻¹=Qᵀ。
assert np.allclose(vectors.T @ vectors, np.eye(2))
# 核对“特征值之和等于迹”；迹在换基后保持不变。
assert np.isclose(values.sum(), np.trace(A))
# 核对“特征值之积等于行列式”；在二维中该乘积也是面积缩放倍数。
assert np.isclose(values.prod(), np.linalg.det(A))
# 用从 0 到 2π 的 200 个角度描出单位圆；含两个端点，曲线自然闭合。
theta = np.linspace(0, 2 * np.pi, 200)
# vstack 把横、纵坐标叠成 (2,200)，A@(2,200) 同时变换 200 个列向量。
# 由于 A 对称且正定，圆变为沿特征向量方向、半轴长度为特征值的椭圆。
circle = np.vstack([np.cos(theta), np.sin(theta)])
ellipse = A @ circle
# 创建一张图及一个坐标轴；figsize 的单位是英寸。
fig, ax = plt.subplots(figsize=(4.5, 4))
# 星号把两行坐标展开为 plot(x,y)；黑色虚线给出变换前的单位圆。
ax.plot(*circle, "k--", label="单位圆")
# 将全部变换后的点连成曲线，与单位圆比较伸缩方向和幅度。
ax.plot(*ellipse, label="A 映射后的圆")
# vectors.T 的每一行恰好是原矩阵的一列，因此可逐对取出 (λ_j,v_j)。
# 这里转置是为了 Python 的逐行迭代习惯，并未改变特征向量按列存储的约定。
for value, vector in zip(values, vectors.T):
    # 把 -λv 和 +λv 作为两列，组成 (2,2) 的线段端点数组，画出完整的椭圆轴。
    endpoints = value * np.column_stack([-vector, vector])
    # 图例显示该方向的伸缩倍数；3 位小数只影响图例文本。
    ax.plot(*endpoints, label=f"主轴 {value:.3f}")
# 横纵坐标使用相同长度单位，避免把圆因画布比例画成椭圆；淡网格辅助读图。
ax.set_aspect("equal")
ax.grid(alpha=0.3)
ax.legend()
# 标题说明本图展示的是对称正定线性变换。
ax.set_title("对称正定矩阵的变换")
# 自动调整边距，以 180 dpi 保存，再关闭图对象以释放内存。
fig.tight_layout()
fig.savefig("e25_eigen.png", dpi=180, bbox_inches="tight")
plt.close(fig)