# 教材位置：第7章，第7.5节；图 7-1、图 7-3。
# 单位、变量与形状：矩阵 A 为 (2,2)，解与右端为 (2,)；本例数值无量纲，条件数衡量相对敏感
# 度。
# 学习目标：求解线性方程组 Ax=b，并用近奇异矩阵观察输入扰动的放大。
# 核心数学关系：可逆方阵给出唯一解；条件数 kappa(A)=||A||_2*||A^{-1}||_2
# 描述最不利方向上的相对误差放大上界，固定 A 时有
# ||delta x||/||x|| <= kappa(A)*||delta b||/||b||（精确算术下）。
# 数据与输出：先解一个 2×2 方程组，再保持病态矩阵 C 不变而轻微改变右端，
# 打印两组解、矩阵性质和相对变化，保存 e23_solve.png，并用 e23_rank.png解释秩与零空间。
# 运行方式：安装 NumPy、Matplotlib 后，运行 python e23_solve.py；
# 图片写入运行命令时的当前工作目录。

import numpy as np
import matplotlib.pyplot as plt
# 在创建第一幅图之前配置可跨系统选择的中文字体。
from plot_style import configure_chinese_plot

configure_chinese_plot()
# A 形状 (2,2)，每一行是一条方程的系数；b 形状 (2,) 是对应右端项。
# 方程为 x+2y=9、3x-y=5，待求解向量包含未知量 x、y。
A = np.array([[1.0, 2.0], [3.0, -1.0]])
b = np.array([9.0, 5.0])
# solve 直接求解 Ax=b，返回形状 (2,)；无需显式计算 inv(A)@b，
# 直接求解通常更高效、数值上也更合适。本接口要求 A 为方阵且非奇异。
x = np.linalg.solve(A, b)
# round 仅控制显示的小数位，x 本身仍保留完整浮点结果供后续检验。
print("方程组的解 =", np.round(x, 6))
# rank 是数值秩，det 是行列式；cond 默认给出 2-范数条件数。
# 满秩说明在精确模型下可唯一求解，条件数则进一步描述解对扰动有多敏感；
# 行列式的绝对大小还受缩放影响，不能单独用它衡量数值稳定性。
print("秩、行列式、条件数 =", np.linalg.matrix_rank(A),
      round(np.linalg.det(A), 6), round(np.linalg.cond(A), 6))
# 同时检查解析解和代回方程的残差；allclose 允许浮点舍入误差。
assert np.allclose(x, [19 / 7, 22 / 7]) and np.allclose(A @ x, b)
# C 的两行几乎相同，仍可逆但接近奇异；两条约束几乎平行，交点很敏感。
# “病态”描述问题本身的敏感性，不能简单归因于求解算法出错。
C = np.array([[1.0, 1.0], [1.0, 1.0001]])
# 只让第二个右端分量增加 0.0001，系数矩阵 C 完全不变。
b1 = np.array([2.0, 2.0])
b2 = np.array([2.0, 2.0001])
# 相同的近奇异矩阵配上略有不同的右端，解分别接近 [2,0] 与 [1,1]。
x1 = np.linalg.solve(C, b1)
x2 = np.linalg.solve(C, b2)
# 使用欧氏范数衡量整体变化，并分别以原右端和原解的范数归一化，
# 得到无单位的相对变化，便于比较不同量纲/规模的问题。
relative_b = np.linalg.norm(b2 - b1) / np.linalg.norm(b1)
relative_x = np.linalg.norm(x2 - x1) / np.linalg.norm(x1)
# 条件数很大意味着存在可被显著放大的扰动方向，但并非每个扰动都取到上界。
print("扰动前后的解 =", x1, x2, "；条件数 =",
      round(np.linalg.cond(C), 2))
print("右端 b、解 x 的相对变化 =", round(relative_b, 8), round(relative_x, 6))
# 检查两次解确实相距明显；这与右端变化很小并不矛盾。
assert np.allclose(x1, [2, 0]) and np.allclose(x2, [1, 1])
# 对“仅扰动 b”检查条件数上界，额外的 1e-10 容差吸收浮点舍入影响。
# 上界是最坏情形的保证，不要求相对解变化恰好等于该上界。
assert relative_x <= np.linalg.cond(C) * relative_b + 1e-10
# 为秩与零空间的概念配套作图：满秩矩阵直接复用本例开头的 A。
# rank_deficient 是另一个明确命名的秩亏矩阵；它与病态但可逆的 C 不是同一对象。
rank_deficient = np.array([[1.0, 1.0], [0.0, 0.0]])
null_vector = np.array([1.0, -1.0])
assert np.linalg.matrix_rank(A) == 2
assert np.linalg.matrix_rank(rank_deficient) == 1
assert np.allclose(rank_deficient @ null_vector, 0)
# 200个单位圆列向量由同一矩阵统一映射；满秩保留面积，秩亏只留下一个方向。
angles = np.linspace(0, 2 * np.pi, 200)
unit_circle = np.vstack([np.cos(angles), np.sin(angles)])
fig, axes = plt.subplots(1, 2, figsize=(6.3, 3.0))
for ax, matrix, title in zip(axes, [A, rank_deficient],
                           ["原矩阵 A：满秩", "秩亏矩阵：秩为 1"]):
    ax.plot(*unit_circle, "k--", alpha=0.5, label="单位圆")
    ax.plot(*(matrix @ unit_circle), lw=1.8, label="映射后的圆")
    ax.axhline(0, color="gray", lw=0.5)
    ax.axvline(0, color="gray", lw=0.5)
    ax.set_aspect("equal")
    ax.grid(alpha=0.25)
    ax.set_title(title, fontsize=10)
    ax.set(xlim=(-3.6, 3.6), ylim=(-3.6, 3.6), xlabel="$x_1$", ylabel="$x_2$")
    ax.legend(fontsize=7.8)
# z=(1,-1)落在rank_deficient的零空间；图中原点就是它变换后的像。
axes[1].annotate("", xy=null_vector, xytext=(0, 0),
                 arrowprops={"arrowstyle": "->", "color": "tab:green"})
axes[1].text(1.05, -1.3, "$z=(1,-1)^T$", fontsize=8, color="tab:green")
axes[1].scatter([0], [0], color="tab:orange", s=22, zorder=4)
axes[1].text(0.03, 0.04, "秩=1，零化度=1",
             transform=axes[1].transAxes, fontsize=8)
fig.tight_layout()
fig.savefig("e23_rank.png", dpi=180, bbox_inches="tight")
plt.close(fig)

# 生成横坐标网格，用二维几何图展示两个很接近的约束方向。
xx = np.linspace(-1, 3, 200)
fig, axes = plt.subplots(1, 2, figsize=(6.3, 3.1))
# 第一条方程 x+y=2 是两个系统共用的约束。
# 另外同时画出原系统与扰动系统的第二条约束，避免把原解误当成另外两线的交点。
def constraint_lines(horizontal):
    # 每条线均从实际 C 与 b1/b2 求出；没有重新输入另一套近奇异系数。
    common = (b1[0] - C[0, 0] * horizontal) / C[0, 1]
    original = (b1[1] - C[1, 0] * horizontal) / C[1, 1]
    perturbed = (b2[1] - C[1, 0] * horizontal) / C[1, 1]
    return common, original, perturbed
labels = ["共同约束：$x+y=2$", "原约束：$x+1.0001y=2$",
          "扰动约束：$x+1.0001y=2.0001$"]
colors = ["tab:blue", "tab:green", "tab:orange"]
styles = ["-", ":", "--"]
ax = axes[0]
for values, label, color, style in zip(constraint_lines(xx), labels, colors, styles):
    ax.plot(xx, values, linestyle=style, color=color, label=label)
# 两个点分别代表两套系统的解，zorder=3 使点显示在曲线之上。
ax.scatter([x1[0], x2[0]], [x1[1], x2[1]], color=["tab:green", "tab:orange"], zorder=3)
ax.annotate("原解 (2,0)", xy=x1, xytext=(1.5, 0.65), fontsize=8, color="tab:green")
ax.annotate("扰动解 (1,1)", xy=x2, xytext=(-0.75, 1.45), fontsize=8, color="tab:orange")
# 标注坐标、图例和网格，帮助把代数解读作两条直线的交点。
ax.set(xlabel="$x$", ylabel="$y$", xlim=(-1, 3), ylim=(-1, 3),
    title="约束相近，解相距较远")
ax.legend(fontsize=7.8)
ax.grid(alpha=0.3)
# 围绕(1.5,0.5)按相同倍数放大横纵坐标，三条几乎重合的约束才能看清。
# 这是显示坐标的变换；求解所用 C、b1、b2 均保持原样。
zoom_factor = 1e4
zoom_center = np.array([1.5, 0.5])
zoom_x = np.linspace(-0.6, 0.6, 150)
local_x = zoom_center[0] + zoom_x / zoom_factor
ax = axes[1]
for values, label, color, style in zip(constraint_lines(local_x), labels, colors,
    styles):
    ax.plot(zoom_x, (values - zoom_center[1]) * zoom_factor,
            linestyle=style, color=color, label=label.split(":")[0])
ax.set(xlabel="$10^4(x-1.5)$", ylabel="$10^4(y-0.5)$", title="局部坐标放大 $10^4$")
ax.legend(fontsize=8)
ax.grid(alpha=0.3)
# 调整布局、保存图片并关闭图，便于批量运行时控制资源占用。
fig.tight_layout()
fig.savefig("e23_solve.png", dpi=180, bbox_inches="tight")
plt.close(fig)