第15章：数值线性代数基础

数值计算中的误差与稳定性

浮点数表示与舍入误差

计算机中的实数表示存在精度限制，导致舍入误差的累积。

条件数

矩阵的条件数衡量了线性系统对输入扰动的敏感性： $\kappa(A) = \|A\| \|A^{-1}\|$

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import norm, solve, lu, qr, svd
from scipy.sparse import diags
import time
import warnings
warnings.filterwarnings('ignore')

def numerical_stability_demo():
    """数值稳定性和误差分析演示"""

    print("数值线性代数：稳定性和误差分析")
    print("=" * 60)

    # 1. 机器精度演示
    print("1. 机器精度分析")
    print("-" * 30)

    machine_epsilon = np.finfo(float).eps
    print(f"机器精度 ε: {machine_epsilon:.2e}")

    # 演示浮点数精度问题
    a = 0.1 + 0.1 + 0.1
    b = 0.3
    print(f"0.1 + 0.1 + 0.1 = {a:.17f}")
    print(f"0.3 = {b:.17f}")
    print(f"相等? {a == b}")
    print(f"差值: {abs(a - b):.2e}")

    # 2. 病态矩阵示例
    print(f"\n2. 病态矩阵分析")
    print("-" * 30)

    # Hilbert矩阵（著名的病态矩阵）
    def hilbert_matrix(n):
        """生成n阶Hilbert矩阵"""
        H = np.zeros((n, n))
        for i in range(n):
            for j in range(n):
                H[i, j] = 1.0 / (i + j + 1)
        return H

    # 不同大小的Hilbert矩阵
    sizes = [3, 5, 8, 10]
    condition_numbers = []

    print("Hilbert矩阵条件数分析:")
    for n in sizes:
        H = hilbert_matrix(n)
        cond_num = np.linalg.cond(H)
        condition_numbers.append(cond_num)
        print(f"H_{n}: 条件数 = {cond_num:.2e}")

    # 3. 条件数对求解精度的影响
    print(f"\n3. 条件数对求解精度的影响")
    print("-" * 40)

    n = 6
    H = hilbert_matrix(n)
    x_exact = np.ones(n)  # 精确解
    b_exact = H @ x_exact  # 精确右端项

    # 添加小扰动
    perturbation_levels = [1e-12, 1e-10, 1e-8, 1e-6]
    results = []

    for eps in perturbation_levels:
        # 扰动右端项
        b_perturbed = b_exact + eps * np.random.randn(n)

        # 求解扰动系统
        try:
            x_computed = solve(H, b_perturbed)
            relative_error = norm(x_computed - x_exact) / norm(x_exact)
            results.append((eps, relative_error))

            print(f"扰动水平 {eps:.0e}: 相对误差 = {relative_error:.2e}")
        except:
            print(f"扰动水平 {eps:.0e}: 求解失败")

    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))

    # 条件数增长
    ax1 = axes[0, 0]
    ax1.semilogy(sizes, condition_numbers, 'o-', linewidth=2, markersize=8)
    ax1.set_xlabel('矩阵大小')
    ax1.set_ylabel('条件数')
    ax1.set_title('Hilbert矩阵条件数增长')
    ax1.grid(True, alpha=0.3)

    # Hilbert矩阵可视化
    ax2 = axes[0, 1]
    H_display = hilbert_matrix(8)
    im = ax2.imshow(H_display, cmap='Blues')
    ax2.set_title('8×8 Hilbert矩阵')
    plt.colorbar(im, ax=ax2)

    # 误差放大效应
    ax3 = axes[0, 2]
    if results:
        perturbations, errors = zip(*results)
        ax3.loglog(perturbations, errors, 'ro-', linewidth=2, markersize=8, label='实际误差')

        # 理论误差界
        cond_H = np.linalg.cond(H)
        theoretical_errors = [cond_H * eps / norm(b_exact) for eps in perturbations]
        ax3.loglog(perturbations, theoretical_errors, 'b--', linewidth=2, label='理论误差界')

        ax3.set_xlabel('输入扰动水平')
        ax3.set_ylabel('输出相对误差')
        ax3.set_title('误差放大效应')
        ax3.legend()
        ax3.grid(True, alpha=0.3)

    # 不同矩阵类型的条件数比较
    ax4 = axes[1, 0]

    matrix_types = {
        '单位矩阵': np.eye(5),
        '对角优势': np.diag([10, 5, 3, 2, 1]) + 0.1*np.random.randn(5, 5),
        'Hilbert-5': hilbert_matrix(5),
        '随机矩阵': np.random.randn(5, 5)
    }

    names = []
    conds = []

    for name, A in matrix_types.items():
        if name != 'Hilbert-5':
            A = A @ A.T  # 确保正定
        try:
            cond_A = np.linalg.cond(A)
            names.append(name)
            conds.append(cond_A)
        except:
            pass

    ax4.bar(range(len(names)), conds, alpha=0.7)
    ax4.set_yscale('log')
    ax4.set_xlabel('矩阵类型')
    ax4.set_ylabel('条件数（对数尺度）')
    ax4.set_title('不同矩阵类型的条件数')
    ax4.set_xticks(range(len(names)))
    ax4.set_xticklabels(names, rotation=45)
    ax4.grid(True, alpha=0.3)

    # 舍入误差累积演示
    ax5 = axes[1, 1]

    def sum_series(n, method='forward'):
        """计算级数和的不同方法"""
        if method == 'forward':
            # 正向累加
            result = 0.0
            for i in range(1, n+1):
                result += 1.0 / i**2
        else:  # backward
            # 反向累加
            result = 0.0
            for i in range(n, 0, -1):
                result += 1.0 / i**2
        return result

    n_values = np.logspace(2, 6, 20).astype(int)
    forward_results = []
    backward_results = []

    for n in n_values:
        forward_results.append(sum_series(n, 'forward'))
        backward_results.append(sum_series(n, 'backward'))

    ax5.semilogx(n_values, forward_results, 'r-', label='正向累加')
    ax5.semilogx(n_values, backward_results, 'b-', label='反向累加')
    ax5.axhline(y=np.pi**2/6, color='black', linestyle='--', label='理论值 π²/6')

    ax5.set_xlabel('项数')
    ax5.set_ylabel('级数和')
    ax5.set_title('舍入误差累积：Σ(1/k²)')
    ax5.legend()
    ax5.grid(True, alpha=0.3)

    # 数值稳定性比较
    ax6 = axes[1, 2]

    # 比较不同算法的数值稳定性
    def unstable_algorithm(x):
        """数值不稳定算法示例"""
        return (x**2 - 1) / (x - 1)

    def stable_algorithm(x):
        """数值稳定算法"""
        return x + 1

    x_values = np.linspace(0.99, 1.01, 100)
    x_values = x_values[x_values != 1.0]  # 避免除零

    unstable_results = []
    stable_results = []

    for x in x_values:
        try:
            unstable_results.append(unstable_algorithm(x))
            stable_results.append(stable_algorithm(x))
        except:
            unstable_results.append(np.nan)
            stable_results.append(np.nan)

    ax6.plot(x_values, unstable_results, 'r-', label='不稳定: (x²-1)/(x-1)', linewidth=2)
    ax6.plot(x_values, stable_results, 'b-', label='稳定: x+1', linewidth=2)

    ax6.set_xlabel('x')
    ax6.set_ylabel('函数值')
    ax6.set_title('数值稳定性比较')
    ax6.legend()
    ax6.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

    return condition_numbers

def matrix_factorizations_numerical():
    """矩阵分解的数值实现"""

    print(f"\n矩阵分解的数值方法")
    print("=" * 40)

    # 生成测试矩阵
    np.random.seed(42)
    n = 6
    A_random = np.random.randn(n, n)
    A_spd = A_random @ A_random.T  # 对称正定矩阵
    A_general = np.random.randn(n, n) + 0.1 * np.eye(n)  # 一般矩阵

    print("测试矩阵信息:")
    print(f"随机对称正定矩阵条件数: {np.linalg.cond(A_spd):.2e}")
    print(f"一般矩阵条件数: {np.linalg.cond(A_general):.2e}")

    # 1. LU分解
    print(f"\n1. LU分解")
    print("-" * 20)

    P, L, U = lu(A_general)

    print("LU分解验证:")
    print(f"||PA - LU||: {norm(P @ A_general - L @ U):.2e}")

    # 检查L和U的性质
    print(f"L是下三角: {np.allclose(L, np.tril(L))}")
    print(f"U是上三角: {np.allclose(U, np.triu(U))}")
    print(f"P是置换矩阵: {np.allclose(P @ P.T, np.eye(n))}")

    # 2. QR分解
    print(f"\n2. QR分解")
    print("-" * 20)

    Q, R = qr(A_general)

    print("QR分解验证:")
    print(f"||A - QR||: {norm(A_general - Q @ R):.2e}")
    print(f"Q是正交矩阵: {np.allclose(Q @ Q.T, np.eye(n))}")
    print(f"R是上三角: {np.allclose(R, np.triu(R))}")

    # 3. SVD分解
    print(f"\n3. 奇异值分解")
    print("-" * 20)

    U_svd, s, Vt = svd(A_general)

    print("SVD分解验证:")
    print(f"||A - UΣVᵀ||: {norm(A_general - U_svd @ np.diag(s) @ Vt):.2e}")
    print(f"奇异值: {s}")
    print(f"条件数 (σ_max/σ_min): {s[0]/s[-1]:.2e}")

    # 4. Cholesky分解（仅对正定矩阵）
    print(f"\n4. Cholesky分解")
    print("-" * 20)

    try:
        from scipy.linalg import cholesky
        L_chol = cholesky(A_spd, lower=True)

        print("Cholesky分解验证:")
        print(f"||A - LLᵀ||: {norm(A_spd - L_chol @ L_chol.T):.2e}")
        print(f"L是下三角: {np.allclose(L_chol, np.tril(L_chol))}")

        # Cholesky vs LU计算复杂度比较
        print(f"Cholesky分解节省约50%计算量（对正定矩阵）")

    except Exception as e:
        print(f"Cholesky分解失败: {e}")

    # 可视化矩阵分解
    fig, axes = plt.subplots(3, 4, figsize=(16, 12))

    # 原始矩阵
    ax = axes[0, 0]
    im = ax.imshow(A_general, cmap='RdBu')
    ax.set_title('原矩阵 A')
    plt.colorbar(im, ax=ax)

    # LU分解
    ax = axes[0, 1]
    im = ax.imshow(L, cmap='Blues')
    ax.set_title('L (下三角)')
    plt.colorbar(im, ax=ax)

    ax = axes[0, 2]
    im = ax.imshow(U, cmap='Reds')
    ax.set_title('U (上三角)')
    plt.colorbar(im, ax=ax)

    ax = axes[0, 3]
    im = ax.imshow(P, cmap='Greys')
    ax.set_title('P (置换)')
    plt.colorbar(im, ax=ax)

    # QR分解
    ax = axes[1, 0]
    im = ax.imshow(Q, cmap='RdBu')
    ax.set_title('Q (正交)')
    plt.colorbar(im, ax=ax)

    ax = axes[1, 1]
    im = ax.imshow(R, cmap='Reds')
    ax.set_title('R (上三角)')
    plt.colorbar(im, ax=ax)

    # SVD分解
    ax = axes[1, 2]
    im = ax.imshow(U_svd, cmap='Blues')
    ax.set_title('U (左奇异向量)')
    plt.colorbar(im, ax=ax)

    ax = axes[1, 3]
    ax.bar(range(len(s)), s)
    ax.set_title('奇异值')
    ax.set_xlabel('索引')
    ax.set_ylabel('值')
    ax.grid(True, alpha=0.3)

    # 分解精度比较
    ax = axes[2, 0]
    methods = ['LU', 'QR', 'SVD']
    errors = [
        norm(P @ A_general - L @ U),
        norm(A_general - Q @ R),
        norm(A_general - U_svd @ np.diag(s) @ Vt)
    ]

    bars = ax.bar(methods, errors)
    ax.set_yscale('log')
    ax.set_ylabel('重构误差')
    ax.set_title('分解精度比较')
    ax.grid(True, alpha=0.3)

    # 计算时间比较
    ax = axes[2, 1]

    def time_factorization(A, method, n_trials=10):
        """测试分解方法的计算时间"""
        times = []
        for _ in range(n_trials):
            start_time = time.time()
            if method == 'LU':
                lu(A)
            elif method == 'QR':
                qr(A)
            elif method == 'SVD':
                svd(A)
            end_time = time.time()
            times.append(end_time - start_time)
        return np.mean(times)

    time_results = {}
    for method in ['LU', 'QR', 'SVD']:
        time_results[method] = time_factorization(A_general, method)

    methods = list(time_results.keys())
    times = list(time_results.values())

    ax.bar(methods, times, alpha=0.7)
    ax.set_ylabel('平均时间 (秒)')
    ax.set_title('计算时间比较')
    ax.grid(True, alpha=0.3)

    # 条件数与分解稳定性
    ax = axes[2, 2]

    # 生成不同条件数的矩阵
    condition_numbers = [1e2, 1e4, 1e6, 1e8]
    lu_errors = []
    qr_errors = []

    for cond in condition_numbers:
        # 构造指定条件数的矩阵
        U_temp, _, Vt_temp = svd(np.random.randn(5, 5))
        singular_vals = np.logspace(0, -np.log10(cond), 5)
        A_temp = U_temp @ np.diag(singular_vals) @ Vt_temp

        # 测试分解精度
        P_temp, L_temp, U_temp = lu(A_temp)
        Q_temp, R_temp = qr(A_temp)

        lu_error = norm(P_temp @ A_temp - L_temp @ U_temp) / norm(A_temp)
        qr_error = norm(A_temp - Q_temp @ R_temp) / norm(A_temp)

        lu_errors.append(lu_error)
        qr_errors.append(qr_error)

    ax.loglog(condition_numbers, lu_errors, 'o-', label='LU分解', linewidth=2)
    ax.loglog(condition_numbers, qr_errors, 's-', label='QR分解', linewidth=2)

    ax.set_xlabel('条件数')
    ax.set_ylabel('相对误差')
    ax.set_title('分解稳定性 vs 条件数')
    ax.legend()
    ax.grid(True, alpha=0.3)

    # 空白子图用于说明
    ax = axes[2, 3]
    ax.text(0.1, 0.8, '分解方法特点:', fontsize=12, weight='bold', transform=ax.transAxes)
    ax.text(0.1, 0.7, '• LU: 一般矩阵，O(n³)', fontsize=10, transform=ax.transAxes)
    ax.text(0.1, 0.6, '• QR: 数值稳定，正交化', fontsize=10, transform=ax.transAxes)
    ax.text(0.1, 0.5, '• SVD: 最稳定，最昂贵', fontsize=10, transform=ax.transAxes)
    ax.text(0.1, 0.4, '• Cholesky: 正定矩阵专用', fontsize=10, transform=ax.transAxes)
    ax.text(0.1, 0.2, '选择原则:', fontsize=12, weight='bold', transform=ax.transAxes)
    ax.text(0.1, 0.1, '• 速度: LU > QR > SVD', fontsize=10, transform=ax.transAxes)
    ax.text(0.1, 0.0, '• 稳定性: SVD > QR > LU', fontsize=10, transform=ax.transAxes)
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)
    ax.axis('off')

    plt.tight_layout()
    plt.show()

def iterative_methods():
    """迭代法求解线性方程组"""

    print(f"\n迭代法求解线性方程组")
    print("=" * 40)

    # 生成对角优势矩阵（保证收敛）
    n = 8
    A = np.random.randn(n, n)
    A = A + (n + 1) * np.eye(n)  # 使其对角占优

    x_true = np.random.randn(n)
    b = A @ x_true

    print(f"系统规模: {n} × {n}")
    print(f"矩阵条件数: {np.linalg.cond(A):.2e}")

    # 直接法求解（作为参考）
    x_direct = solve(A, b)
    print(f"直接法误差: {norm(x_direct - x_true):.2e}")

    def jacobi_method(A, b, x0, max_iter=100, tol=1e-10):
        """雅可比迭代法"""
        n = len(b)
        x = x0.copy()
        D = np.diag(np.diag(A))
        R = A - D

        residuals = []
        errors = []

        for k in range(max_iter):
            x_new = solve(D, b - R @ x)

            residual = norm(A @ x_new - b)
            error = norm(x_new - x_true)
            residuals.append(residual)
            errors.append(error)

            if residual < tol:
                return x_new, k+1, residuals, errors

            x = x_new

        return x, max_iter, residuals, errors

    def gauss_seidel_method(A, b, x0, max_iter=100, tol=1e-10):
        """高斯-赛德尔迭代法"""
        n = len(b)
        x = x0.copy()

        residuals = []
        errors = []

        for k in range(max_iter):
            x_new = x.copy()

            for i in range(n):
                sum1 = np.dot(A[i, :i], x_new[:i])
                sum2 = np.dot(A[i, i+1:], x[i+1:])
                x_new[i] = (b[i] - sum1 - sum2) / A[i, i]

            residual = norm(A @ x_new - b)
            error = norm(x_new - x_true)
            residuals.append(residual)
            errors.append(error)

            if residual < tol:
                return x_new, k+1, residuals, errors

            x = x_new

        return x, max_iter, residuals, errors

    def sor_method(A, b, x0, omega=1.2, max_iter=100, tol=1e-10):
        """超松弛迭代法"""
        n = len(b)
        x = x0.copy()

        residuals = []
        errors = []

        for k in range(max_iter):
            x_new = x.copy()

            for i in range(n):
                sum1 = np.dot(A[i, :i], x_new[:i])
                sum2 = np.dot(A[i, i+1:], x[i+1:])
                x_gs = (b[i] - sum1 - sum2) / A[i, i]
                x_new[i] = (1 - omega) * x[i] + omega * x_gs

            residual = norm(A @ x_new - b)
            error = norm(x_new - x_true)
            residuals.append(residual)
            errors.append(error)

            if residual < tol:
                return x_new, k+1, residuals, errors

            x = x_new

        return x, max_iter, residuals, errors

    # 运行各种迭代方法
    x0 = np.zeros(n)  # 初始猜测

    print(f"\n迭代方法比较:")
    print("-" * 30)

    methods = {
        'Jacobi': jacobi_method,
        'Gauss-Seidel': gauss_seidel_method,
        'SOR (ω=1.2)': lambda A, b, x0, **kwargs: sor_method(A, b, x0, omega=1.2, **kwargs)
    }

    results = {}

    for name, method in methods.items():
        x_sol, iterations, residuals, errors = method(A, b, x0)
        final_error = norm(x_sol - x_true)

        results[name] = {
            'solution': x_sol,
            'iterations': iterations,
            'residuals': residuals,
            'errors': errors,
            'final_error': final_error
        }

        print(f"{name:12}: {iterations:3d} 次迭代，最终误差 = {final_error:.2e}")

    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))

    # 收敛历史
    ax1 = axes[0, 0]
    for name, result in results.items():
        iterations_list = range(1, len(result['residuals']) + 1)
        ax1.semilogy(iterations_list, result['residuals'], 'o-',
                    linewidth=2, markersize=4, label=name)

    ax1.set_xlabel('迭代次数')
    ax1.set_ylabel('残差 ||Ax - b||')
    ax1.set_title('残差收敛历史')
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # 误差收敛历史
    ax2 = axes[0, 1]
    for name, result in results.items():
        iterations_list = range(1, len(result['errors']) + 1)
        ax2.semilogy(iterations_list, result['errors'], 'o-',
                    linewidth=2, markersize=4, label=name)

    ax2.set_xlabel('迭代次数')
    ax2.set_ylabel('误差 ||x - x_true||')
    ax2.set_title('误差收敛历史')
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    # 收敛速度比较
    ax3 = axes[0, 2]
    method_names = list(results.keys())
    iterations_count = [results[name]['iterations'] for name in method_names]
    final_errors = [results[name]['final_error'] for name in method_names]

    bars = ax3.bar(range(len(method_names)), iterations_count, alpha=0.7)
    ax3.set_xlabel('方法')
    ax3.set_ylabel('收敛所需迭代次数')
    ax3.set_title('收敛速度比较')
    ax3.set_xticks(range(len(method_names)))
    ax3.set_xticklabels(method_names, rotation=45)
    ax3.grid(True, alpha=0.3)

    # 添加数值标签
    for bar, count in zip(bars, iterations_count):
        ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.5,
                f'{count}', ha='center', va='bottom')

    # SOR参数优化
    ax4 = axes[1, 0]
    omega_values = np.linspace(0.5, 1.95, 30)
    sor_iterations = []

    for omega in omega_values:
        try:
            _, iters, _, _ = sor_method(A, b, x0, omega=omega, max_iter=200)
            sor_iterations.append(iters)
        except:
            sor_iterations.append(200)  # 未收敛

    ax4.plot(omega_values, sor_iterations, 'b-', linewidth=2)
    optimal_idx = np.argmin(sor_iterations)
    optimal_omega = omega_values[optimal_idx]
    ax4.plot(optimal_omega, sor_iterations[optimal_idx], 'ro', markersize=10,
             label=f'最优 ω = {optimal_omega:.2f}')

    ax4.set_xlabel('松弛因子 ω')
    ax4.set_ylabel('收敛迭代次数')
    ax4.set_title('SOR方法参数优化')
    ax4.legend()
    ax4.grid(True, alpha=0.3)

    # 谱半径分析
    ax5 = axes[1, 1]

    # 计算迭代矩阵的谱半径
    D = np.diag(np.diag(A))
    L = np.tril(A, k=-1)
    U = np.triu(A, k=1)

    # Jacobi迭代矩阵
    T_jacobi = -solve(D, L + U)
    rho_jacobi = max(abs(np.linalg.eigvals(T_jacobi)))

    # Gauss-Seidel迭代矩阵
    T_gs = -solve(D + L, U)
    rho_gs = max(abs(np.linalg.eigvals(T_gs)))

    methods_spectral = ['Jacobi', 'Gauss-Seidel']
    spectral_radii = [rho_jacobi, rho_gs]

    bars = ax5.bar(methods_spectral, spectral_radii, alpha=0.7)
    ax5.axhline(y=1, color='red', linestyle='--', label='收敛界限')
    ax5.set_ylabel('谱半径')
    ax5.set_title('迭代矩阵谱半径')
    ax5.legend()
    ax5.grid(True, alpha=0.3)

    # 添加数值标签
    for bar, rho in zip(bars, spectral_radii):
        ax5.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01,
                f'{rho:.3f}', ha='center', va='bottom')

    # 解的收敛轨迹（选择前两个分量）
    ax6 = axes[1, 2]

    for name, result in results.items():
        if name == 'Jacobi':  # 只显示一种方法的轨迹，避免图形混乱
            # 重新运行以获取所有中间解
            x = x0.copy()
            trajectory_x = [x[0]]
            trajectory_y = [x[1]]

            D = np.diag(np.diag(A))
            R = A - D

            for k in range(min(50, result['iterations'])):
                x = solve(D, b - R @ x)
                trajectory_x.append(x[0])
                trajectory_y.append(x[1])

            ax6.plot(trajectory_x, trajectory_y, 'b-o', markersize=4,
                    alpha=0.7, label='Jacobi轨迹')

    ax6.plot(x_true[0], x_true[1], 'rs', markersize=10, label='真解')
    ax6.plot(x0[0], x0[1], 'go', markersize=8, label='初始值')

    ax6.set_xlabel('x₁')
    ax6.set_ylabel('x₂')
    ax6.set_title('解的收敛轨迹')
    ax6.legend()
    ax6.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

    return results

def eigenvalue_algorithms():
    """特征值计算的数值方法"""

    print(f"\n特征值计算的数值方法")
    print("=" * 40)

    # 生成测试矩阵
    np.random.seed(42)
    n = 6
    A_random = np.random.randn(n, n)
    A_symmetric = A_random + A_random.T  # 对称矩阵
    A_general = A_random + 0.1 * np.eye(n)  # 一般矩阵

    print("测试矩阵:")
    print(f"对称矩阵 A_sym 的条件数: {np.linalg.cond(A_symmetric):.2e}")
    print(f"一般矩阵 A_gen 的条件数: {np.linalg.cond(A_general):.2e}")

    # 真实特征值（用作参考）
    true_eigenvals_sym = np.linalg.eigvals(A_symmetric)
    true_eigenvals_gen = np.linalg.eigvals(A_general)

    # 按模长排序
    true_eigenvals_sym = true_eigenvals_sym[np.argsort(np.abs(true_eigenvals_sym))[::-1]]
    true_eigenvals_gen = true_eigenvals_gen[np.argsort(np.abs(true_eigenvals_gen))[::-1]]

    print(f"\n对称矩阵真实特征值: {true_eigenvals_sym}")
    print(f"一般矩阵真实特征值: {true_eigenvals_gen}")

    def power_method(A, max_iter=100, tol=1e-10):
        """幂法求主特征值"""
        n = A.shape[0]
        x = np.random.randn(n)
        x = x / norm(x)

        eigenval_history = []
        eigenvec_history = []

        for k in range(max_iter):
            y = A @ x
            eigenval = x.T @ y
            x_new = y / norm(y)

            eigenval_history.append(eigenval)
            eigenvec_history.append(x_new.copy())

            if k > 0 and abs(eigenval_history[-1] - eigenval_history[-2]) < tol:
                return eigenval, x_new, k+1, eigenval_history

            x = x_new

        return eigenval_history[-1], x, max_iter, eigenval_history

    def inverse_power_method(A, sigma, max_iter=100, tol=1e-10):
        """反幂法求接近σ的特征值"""
        n = A.shape[0]
        x = np.random.randn(n)
        x = x / norm(x)

        A_shifted = A - sigma * np.eye(n)

        eigenval_history = []

        for k in range(max_iter):
            try:
                y = solve(A_shifted, x)
                mu = x.T @ y
                x_new = y / norm(y)

                eigenval = 1/mu + sigma  # 逆变换回原特征值

                eigenval_history.append(eigenval)

                if k > 0 and abs(eigenval_history[-1] - eigenval_history[-2]) < tol:
                    return eigenval, x_new, k+1, eigenval_history

                x = x_new
            except:
                break

        return eigenval_history[-1] if eigenval_history else sigma, x, k+1, eigenval_history

    def qr_algorithm_basic(A, max_iter=100, tol=1e-10):
        """基础QR算法"""
        A_k = A.copy()
        eigenval_history = []

        for k in range(max_iter):
            Q, R = qr(A_k)
            A_k = R @ Q

            # 记录对角元素（特征值的近似）
            eigenvals = np.diag(A_k)
            eigenval_history.append(eigenvals.copy())

            # 检查收敛（下三角部分接近零）
            if k > 5:
                lower_triangular_norm = norm(np.tril(A_k, k=-1))
                if lower_triangular_norm < tol:
                    break

        return np.diag(A_k), A_k, k+1, eigenval_history

    # 运行各种方法
    print(f"\n幂法求主特征值:")
    print("-" * 25)

    # 对称矩阵
    eigenval_pm_sym, eigenvec_pm_sym, iter_pm_sym, hist_pm_sym = power_method(A_symmetric)
    error_pm_sym = abs(eigenval_pm_sym - true_eigenvals_sym[0])

    print(f"对称矩阵:")
    print(f"  幂法结果: {eigenval_pm_sym:.6f} ({iter_pm_sym} 次迭代)")
    print(f"  真实值: {true_eigenvals_sym[0]:.6f}")
    print(f"  误差: {error_pm_sym:.2e}")

    # 一般矩阵
    eigenval_pm_gen, eigenvec_pm_gen, iter_pm_gen, hist_pm_gen = power_method(A_general)
    error_pm_gen = abs(eigenval_pm_gen - true_eigenvals_gen[0])

    print(f"一般矩阵:")
    print(f"  幂法结果: {eigenval_pm_gen:.6f} ({iter_pm_gen} 次迭代)")
    print(f"  真实值: {true_eigenvals_gen[0]:.6f}")
    print(f"  误差: {error_pm_gen:.2e}")

    # 反幂法
    print(f"\n反幂法求特定特征值:")
    print("-" * 30)

    # 选择一个目标值
    target = 0.5
    eigenval_ipm, eigenvec_ipm, iter_ipm, hist_ipm = inverse_power_method(A_symmetric, target)

    # 找到最接近的真实特征值
    closest_true = true_eigenvals_sym[np.argmin(np.abs(true_eigenvals_sym - target))]
    error_ipm = abs(eigenval_ipm - closest_true)

    print(f"目标值: {target}")
    print(f"反幂法结果: {eigenval_ipm:.6f} ({iter_imp} 次迭代)" if 'iter_ipm' in locals() else f"反幂法结果: {eigenval_ipm:.6f}")
    print(f"最接近的真实特征值: {closest_true:.6f}")
    print(f"误差: {error_ipm:.2e}")

    # QR算法
    print(f"\nQR算法求所有特征值:")
    print("-" * 30)

    eigenvals_qr_sym, A_final_sym, iter_qr_sym, hist_qr_sym = qr_algorithm_basic(A_symmetric)
    eigenvals_qr_gen, A_final_gen, iter_qr_gen, hist_qr_gen = qr_algorithm_basic(A_general)

    print(f"对称矩阵 ({iter_qr_sym} 次迭代):")
    print(f"  QR结果: {np.sort(eigenvals_qr_sym)[::-1]}")
    print(f"  真实值: {np.sort(true_eigenvals_sym)[::-1]}")
    print(f"  误差: {norm(np.sort(eigenvals_qr_sym) - np.sort(true_eigenvals_sym)):.2e}")

    print(f"一般矩阵 ({iter_qr_gen} 次迭代):")
    print(f"  QR结果: {np.sort(np.real(eigenvals_qr_gen))[::-1]}")
    print(f"  真实值: {np.sort(np.real(true_eigenvals_gen))[::-1]}")

    # 可视化
    fig, axes = plt.subplots(3, 3, figsize=(18, 18))

    # 幂法收敛历史
    ax1 = axes[0, 0]
    ax1.semilogy(range(1, len(hist_pm_sym)+1),
                [abs(val - true_eigenvals_sym[0]) for val in hist_pm_sym],
                'b-o', linewidth=2, markersize=4, label='对称矩阵')

    ax1.semilogy(range(1, len(hist_pm_gen)+1),
                [abs(val - true_eigenvals_gen[0]) for val in hist_pm_gen],
                'r-s', linewidth=2, markersize=4, label='一般矩阵')

    ax1.set_xlabel('迭代次数')
    ax1.set_ylabel('误差')
    ax1.set_title('幂法收敛历史')
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # 反幂法收敛历史
    ax2 = axes[0, 1]
    if hist_ipm:
        ax2.semilogy(range(1, len(hist_ipm)+1),
                    [abs(val - closest_true) for val in hist_ipm],
                    'g-o', linewidth=2, markersize=4)

        ax2.set_xlabel('迭代次数')
        ax2.set_ylabel('误差')
        ax2.set_title(f'反幂法收敛历史 (目标={target})')
        ax2.grid(True, alpha=0.3)

    # QR算法收敛历史（显示最大非对角元素）
    ax3 = axes[0, 2]

    def off_diagonal_norm(eigenval_matrix):
        """计算矩阵的非对角元素范数"""
        return norm(eigenval_matrix - np.diag(np.diag(eigenval_matrix)))

    if len(hist_qr_sym) > 1:
        # 这里需要重新运行QR算法来获取矩阵历史
        # 简化显示：只显示特征值的变化
        max_changes_sym = []
        for i in range(1, len(hist_qr_sym)):
            change = max(abs(hist_qr_sym[i] - hist_qr_sym[i-1]))
            max_changes_sym.append(change)

        if max_changes_sym:
            ax3.semilogy(range(1, len(max_changes_sym)+1), max_changes_sym,
                        'b-o', linewidth=2, markersize=4, label='对称矩阵')

    ax3.set_xlabel('迭代次数')
    ax3.set_ylabel('最大特征值变化')
    ax3.set_title('QR算法收敛历史')
    ax3.legend()
    ax3.grid(True, alpha=0.3)

    # 特征值比较
    ax4 = axes[1, 0]
    x_pos = range(len(true_eigenvals_sym))
    width = 0.35

    ax4.bar([x - width/2 for x in x_pos], np.sort(true_eigenvals_sym)[::-1], width,
           label='真实值', alpha=0.7)
    ax4.bar([x + width/2 for x in x_pos], np.sort(eigenvals_qr_sym)[::-1], width,
           label='QR算法', alpha=0.7)

    ax4.set_xlabel('特征值索引')
    ax4.set_ylabel('特征值')
    ax4.set_title('对称矩阵特征值比较')
    ax4.legend()
    ax4.grid(True, alpha=0.3)

    # 计算时间比较
    ax5 = axes[1, 1]

    # 简化的时间测试
    methods = ['Power Method', 'QR Algorithm', 'NumPy (eig)']
    times = []

    # 幂法时间（估算）
    times.append(iter_pm_sym * 1e-4)  # 估算每次迭代时间

    # QR算法时间（估算）
    times.append(iter_qr_sym * 1e-3)  # QR每次迭代更昂贵

    # NumPy时间
    start_time = time.time()
    for _ in range(10):
        np.linalg.eig(A_symmetric)
    numpy_time = (time.time() - start_time) / 10
    times.append(numpy_time)

    bars = ax5.bar(methods, times, alpha=0.7)
    ax5.set_ylabel('时间 (秒)')
    ax5.set_title('计算时间比较 (估算)')
    ax5.set_xticklabels(methods, rotation=45)

    for bar, time_val in zip(bars, times):
        ax5.text(bar.get_x() + bar.get_width()/2, bar.get_height(),
                f'{time_val:.2e}', ha='center', va='bottom')

    # 特征向量可视化（前两个分量）
    ax6 = axes[1, 2]

    # 比较主特征向量
    true_eigenvec = np.linalg.eig(A_symmetric)[1][:, 0]

    # 确保符号一致
    if np.dot(eigenvec_pm_sym, true_eigenvec) < 0:
        eigenvec_pm_sym = -eigenvec_pm_sym

    ax6.arrow(0, 0, true_eigenvec[0], true_eigenvec[1],
             head_width=0.05, head_length=0.05, fc='blue', ec='blue',
             linewidth=3, label='真实特征向量')
    ax6.arrow(0, 0, eigenvec_pm_sym[0], eigenvec_pm_sym[1],
             head_width=0.05, head_length=0.05, fc='red', ec='red',
             linewidth=3, alpha=0.7, label='幂法结果')

    ax6.set_xlim(-1.2, 1.2)
    ax6.set_ylim(-1.2, 1.2)
    ax6.set_xlabel('v₁')
    ax6.set_ylabel('v₂')
    ax6.set_title('主特征向量比较')
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    ax6.set_aspect('equal')

    # 算法复杂度比较
    ax7 = axes[2, 0]

    matrix_sizes = [10, 20, 50, 100]
    complexities = {
        'Power Method': [n for n in matrix_sizes],  # O(n²) per iteration
        'QR Algorithm': [n**3 for n in matrix_sizes],  # O(n³)
        'Direct (eig)': [n**3 for n in matrix_sizes]   # O(n³)
    }

    for method, complexity in complexities.items():
        ax7.loglog(matrix_sizes, complexity, 'o-', linewidth=2, label=method)

    ax7.set_xlabel('矩阵大小 n')
    ax7.set_ylabel('计算复杂度')
    ax7.set_title('算法复杂度比较')
    ax7.legend()
    ax7.grid(True, alpha=0.3)

    # 收敛速度比较
    ax8 = axes[2, 1]

    # 模拟不同特征值间隙对幂法收敛的影响
    gaps = [0.1, 0.3, 0.5, 0.7, 0.9]  # λ₂/λ₁ 比值
    convergence_rates = []

    for gap in gaps:
        # 理论收敛率 ≈ |λ₂/λ₁|
        rate = gap
        convergence_rates.append(rate)

    ax8.plot(gaps, convergence_rates, 'bo-', linewidth=2, markersize=8)
    ax8.set_xlabel('特征值间隙 (λ₂/λ₁)')
    ax8.set_ylabel('收敛率')
    ax8.set_title('幂法收敛率 vs 特征值间隙')
    ax8.grid(True, alpha=0.3)

    # 方法选择指南
    ax9 = axes[2, 2]
    ax9.text(0.1, 0.9, '特征值算法选择指南', fontsize=14, weight='bold', transform=ax9.transAxes)

    guidelines = [
        '单个主特征值：',
        '• 幂法 - 简单、内存效率高',
        '• 适用于大型稀疏矩阵',
        '',
        '特定特征值：',
        '• 反幂法 - 收敛快',
        '• 需要近似值',
        '',
        '所有特征值：',
        '• QR算法 - 数值稳定',
        '• 对称矩阵优化版本',
        '',
        '大规模问题：',
        '• Lanczos/Arnoldi方法',
        '• Krylov子空间方法'
    ]

    for i, text in enumerate(guidelines):
        y_pos = 0.8 - i * 0.05
        if text.startswith('•'):
            ax9.text(0.15, y_pos, text, fontsize=10, transform=ax9.transAxes)
        elif text == '':
            continue
        else:
            ax9.text(0.1, y_pos, text, fontsize=11, weight='bold', transform=ax9.transAxes)

    ax9.set_xlim(0, 1)
    ax9.set_ylim(0, 1)
    ax9.axis('off')

    plt.tight_layout()
    plt.show()

# 运行所有演示
condition_nums = numerical_stability_demo()
matrix_factorizations_numerical()
iter_results = iterative_methods()
eigenvalue_algorithms()

# 最终总结
def numerical_linear_algebra_summary():
    """数值线性代数总结"""

    print(f"\n数值线性代数总结")
    print("=" * 50)

    summary = {
        "误差分析": {
            "舍入误差": "浮点数精度限制",
            "条件数": "系统对输入扰动的敏感性",
            "数值稳定性": "算法对舍入误差的鲁棒性",
            "后向误差": "实际解对应的扰动问题"
        },

        "矩阵分解": {
            "LU分解": "A = PLU，高效求解",
            "QR分解": "A = QR，数值稳定",
            "SVD分解": "A = UΣVᵀ，最稳定最通用",
            "Cholesky": "A = LLᵀ，正定矩阵专用"
        },

        "线性方程组": {
            "直接法": "分解求解，精确但昂贵",
            "Jacobi": "简单并行，收敛慢",
            "Gauss-Seidel": "收敛更快，串行",
            "SOR": "加速收敛，需调参"
        },

        "特征值问题": {
            "幂法": "主特征值，简单高效",
            "反幂法": "特定特征值，快速收敛",
            "QR算法": "所有特征值，工业标准",
            "Krylov方法": "大规模稀疏矩阵"
        }
    }

    for category, methods in summary.items():
        print(f"\n{category}:")
        print("-" * 30)
        for method, description in methods.items():
            print(f"• {method:15}: {description}")

    print(f"\n实践建议:")
    print("-" * 20)
    recommendations = [
        "• 优先考虑数值稳定性而非速度",
        "• 利用矩阵的特殊结构（对称、稀疏等）",
        "• 监控条件数，警惕病态问题",
        "• 选择合适的分解方法",
        "• 使用成熟的数值库（LAPACK, BLAS）",
        "• 验证结果的数值精度"
    ]

    for rec in recommendations:
        print(rec)

numerical_linear_algebra_summary()

# 完成所有章节
print(f"\n" + "="*60)
print("线性代数完整课程学习笔记生成完成！")
print("="*60)
print(f"已生成章节：")
for i in range(8, 16):
    print(f"• 第{i}章：la-{i:03d}.md")

# 最后更新todo列表
from datetime import datetime
print(f"\n生成完成时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")