CE7453 Numerical Algorithms 期末高分超详细攻略 · 章节导航

本系列笔记覆盖 CE7453 期末考试所有核心模块，每章独立 markdown 文件，内容包括详细原理、公式推导、算法流程、例题全解、常见考点与英文关键词。适合零基础、考前冲刺、查漏补缺。

章节索引

Root Finding（方程求根）
- 二分法、牛顿法、割线法、误差与收敛性、典型例题
- 文件：CE7453_期末高分速成.md（含线性方程组）
Linear Systems（线性方程组）
- 高斯消元、LU分解、Jacobi/Gauss-Seidel、梯度下降、例题
- 文件：CE7453_期末高分速成.md
Bezier/B-spline/Interpolation（曲线与插值）
- Bezier曲线、B-spline、de Casteljau与de Boor算法、插值多项式、例题
- 文件：CE7453_03_Bezier_Bspline_Interpolation.md
Numerical Differentiation & Integration（数值微分与积分）
- 差分公式、梯形法、Simpson法、Romberg、Gaussian Quadrature、例题
- 文件：CE7453_04_Numerical_Diff_Integration.md
Least Squares（最小二乘）
- 正规方程、QR分解、非线性最小二乘、GPS应用、例题
- 文件：CE7453_05_Least_Squares.md
Eigenanalysis（特征值分析）
- 幂迭代、QR算法、Rayleigh商、PageRank、例题
- 文件：CE7453_06_Eigenanalysis.md
Fourier Transform（傅里叶变换）
- DFT、FFT、三角插值、DCT、JPEG应用、例题
- 文件：CE7453_07_Fourier_Transform.md
3D Data Registration（3D数据配准）
- ICP算法、SVD分解、点云配准、例题
- 文件：CE7453_08_3D_Data_Registration.md

使用建议

每章内容均为独立 markdown 文件，可按需查阅或打印。
建议先看导航，按考试大纲查漏补缺。
例题部分适合考前手算练习，英文关键词便于查找原版资料。
如需某一算法的代码实现或更详细推导，可随时补充。

祝你期末高分！

CE7453 Numerical Algorithms 期末高分超详细攻略（第一章：Root Finding）

本章内容：方程求根（Root Finding）
适用对象：零基础/考前冲刺/快速查漏补缺
内容特色：详细原理、公式推导、算法流程、例题全解、英文关键词

1.1 基本概念与背景

什么是方程求根？
给定一个连续函数 $f (x)$ ，我们想找到 $x = r$ 使得 $f (r) = 0$ 。这个 $r$ 就叫做 $f (x)$ 的根（root）。
数学表示：
求解方程 $f (x) = 0$ 的解 $x = r$ 。
实际应用：
- 物理学：平衡点、临界值、交叉点计算
- 工程学：结构稳定性分析、电路设计
- 金融学：期权定价、收益率计算
- 计算机图形学：曲线/曲面求交
- 优化问题：目标函数的极值点（导数为零）
求根方法分类：
- 直接法：代数方程的解析解（如一元二次方程公式）
- 迭代法：通过迭代逼近根（本章重点）
- 混合法：结合多种方法的优点

1.2 主要算法与原理

1.2.1 二分法（Bisection Method）

原理：
如果连续函数 $f (x)$ 在区间 $[a, b]$ 上满足 $f (a)$ 和 $f (b)$ 异号（即 $f (a) f (b) < 0$ ），根据中值定理，区间内必存在至少一个根。
算法流程：
1. 检查 $f (a) f (b) < 0$ ，否则不能用二分法。
2. 计算中点 $c = (a + b) /2$ 。
3. 判断 $f (c)$ 的符号：
  - 如果 $f (c) = 0$ ，找到精确根，返回 $c$ 。
  - 如果 $f (a) f (c) < 0$ ，根在 $[a, c]$ ，令 $b = c$ 。
  - 如果 $f (c) f (b) < 0$ ，根在 $[c, b]$ ，令 $a = c$ 。
4. 重复2-3，直到区间足够小（ $∣ b - a ∣ < ϵ$ ）或迭代次数达到上限。
收敛性分析：
- 每次迭代后，区间长度减半： $∣ b_{n + 1} - a_{n + 1} ∣ = \frac{∣ b _{n} - a _{n} ∣}{2}$
- 经过 $n$ 次迭代后，区间长度： $∣ b_{n} - a_{n} ∣ = \frac{∣ b _{0} - a _{0} ∣}{2 ^{n}}$
- 要达到精度 $ϵ$ ，需要迭代次数： $n \geq lo g_{2} \frac{∣ b _{0} - a _{0} ∣}{ϵ}$
误差估计：
- 若 $c_{n}$ 是第 $n$ 次迭代的中点，则 $∣ c_{n} - r ∣ \leq \frac{∣ b _{0} - a _{0} ∣}{2 ^{n + 1}}$ ，其中 $r$ 是真实根。
- 若要求解精确到小数点后 $d$ 位，需满足 $\frac{∣ b _{0} - a _{0} ∣}{2 ^{n}} < 5 \times 1 0^{- (d + 1)}$
- 精度要求：若要求解精确到小数点后 $p$ 位，即误差小于 $0.5 \times 1 0^{- p}$ ，则需要满足： $\frac{∣ b _{0} - a _{0} ∣}{2 ^{n}} < 0.5 \times 1 0^{- p}$
- 所需迭代次数：解上述不等式可得最小迭代次数 $n$ ： $n > lo g_{2} (\frac{∣ b _{0} - a _{0} ∣}{0.5 \times 1 0 ^{- p}})$ 或更方便计算的形式： $n > \frac{lo g _{10} ( ∣ b _{0} - a _{0} ∣ ) - lo g _{10} ( 0.5 \times 1 0 ^{- p} )}{lo g _{10} ( 2 )}$
优缺点：
- 优点：
  - 一定收敛（全局收敛性，global convergence）
  - 实现简单，稳定可靠
  - 不需要导数信息
- 缺点：
  - 收敛速度慢（线性收敛，收敛阶为1）
  - 每次迭代只减少一位有效数字
  - 只适用于区间端点异号的情况
  - 无法处理重根

伪代码：

function Bisection(f, a, b, tol, max_iter)
    if f(a)*f(b) >= 0 then
        return "Error: 区间端点不异号"
    end if
    
    iter = 0
    while (b-a)/2 > tol and iter < max_iter do
        c = (a+b)/2
        if f(c) == 0 then
            return c  // 找到精确根
        else if f(a)*f(c) < 0 then
            b = c
        else
            a = c
        end if
        iter = iter + 1
    end while
    
    return (a+b)/2  // 返回区间中点作为近似根
end function

公式：
$c_{n} = \frac{a _{n} + b _{n}}{2}$
$∣ b_{n} - a_{n} ∣ = \frac{∣ b _{0} - a _{0} ∣}{2 ^{n}}$
英文关键词：bisection method, bracket, interval halving, binary search, convergence, tolerance, error bound

例题1：二分法求解 $f (x) = x^{3} + x - 1$

问题：用二分法求 $f (x) = x^{3} + x - 1$ 在 $[0, 1]$ 内的根，精度要求为 $1 0^{- 2}$ 。

分析： $f (0) = - 1$ , $f (1) = 1$ 。 $f (0) f (1) < 0$ ，符合二分法使用条件。初始区间 $[a_{0}, b_{0}] = [0, 1]$ ，区间长度为 1。要求精度 $ϵ = 1 0^{- 2} = 0.01$ 。需要满足 $\frac{∣ b _{0} - a _{0} ∣}{2 ^{n}} < ϵ$ ，即 $\frac{1}{2 ^{n}} < 0.01$ 。 $2^{n} > 100$ 。因为 $2^{6} = 64$ , $2^{7} = 128$ ，所以至少需要 $n = 7$ 次迭代。

迭代过程：

n	$a_{n}$	$b_{n}$	$c_{n} = (a_{n} + b_{n}) /2$	$f (c_{n})$	区间	长度
0	0	1	0.5	-0.375	[0.5, 1]	0.5
1	0.5	1	0.75	0.171875	[0.5, 0.75]	0.25
2	0.5	0.75	0.625	-0.130859	[0.625, 0.75]	0.125
3	0.625	0.75	0.6875	0.012451	[0.625, 0.6875]	0.0625
4	0.625	0.6875	0.65625	-0.061012	[0.65625, 0.6875]	0.03125
5	0.65625	0.6875	0.671875	-0.024755	[0.671875, 0.6875]	0.015625
6	0.671875	0.6875	0.6796875	-0.006270	[0.6796875, 0.6875]	0.0078125

结果：经过7次迭代，得到近似根 $x \approx c_{6} = 0.6796875$ 。此时区间长度为 $0.0078125 < 0.01$ ，满足精度要求。

Python代码示例：

def is_close(a, b, rtol=1e-5, atol=1e-8):
    """
    检查两个数值是否足够接近，考虑浮点数精度问题
    """
    return np.abs(a - b) <= (atol + rtol * np.abs(b))

def test_example1_bisection():
    """
    验证例题1：用二分法求 f(x)=x^3+x-1 在 [0,1] 内的根，精度 10^-2
    文档中的答案：x ≈ 0.6796875
    """
    f = lambda x: x**3 + x - 1
    
    # 二分法实现
    def bisection_method(f, a, b, tol=1e-10, max_iter=100):
        if f(a) * f(b) >= 0:
            raise ValueError("区间端点函数值必须异号")
        
        iter_count = 0
        while (b - a) / 2 > tol and iter_count < max_iter:
            c = (a + b) / 2
            if f(c) == 0:
                return c
            elif f(a) * f(c) < 0:
                b = c
            else:
                a = c
            iter_count += 1
        
        return (a + b) / 2
    
    # 用自己实现的二分法计算
    root_bisection = bisection_method(f, 0, 1, tol=1e-2)
    # 用scipy的方法计算精确值
    root_exact = optimize.root_scalar(f, bracket=[0, 1], method='brentq').root
    
    # 文档中给出的答案
    root_document = 0.6796875
    
    # 验证自己实现的方法是否与文档答案一致
    # 这里我们只展示计算过程，验证部分在test文件中执行
    print(f"二分法计算结果: {root_bisection}")
    print(f"文档答案: {root_document}")
    print(f"Scipy精确解: {root_exact}")
    print(f"函数在文档答案处的值: {f(root_document)}")

# 运行示例 (需要 import numpy as np 和 from scipy import optimize)
# test_example1_bisection()

1.2.2 牛顿法（Newton's Method）

原理：
利用函数的局部线性近似（泰勒展开）逐步逼近根。在当前点 $x_{n}$ 处作切线，切线与 $x$ 轴的交点作为下一个近似值 $x_{n + 1}$ 。
数学推导：
在 $x_{n}$ 处展开 $f (x)$ 的一阶泰勒级数： $f (x) \approx f (x_{n}) + f^{'} (x_{n}) (x - x_{n})$

令 $f (x) = 0$ ，解得： $x = x_{n} - \frac{f ( x _{n} )}{f ^{'} ( x _{n} )}$

这就是牛顿法的迭代公式： $x_{n + 1} = x_{n} - \frac{f ( x _{n} )}{f ^{'} ( x _{n} )}$
几何解释：
每次迭代相当于用切线近似函数，切线与 $x$ 轴的交点作为下一个近似值。
几何解释：在点 $(x_{n}, f (x_{n}))$ 处作函数 $f (x)$ 的切线，切线的方程为 $y - f (x_{n}) = f^{'} (x_{n}) (x - x_{n})$ 。令 $y = 0$ ，解出切线与 $x$ 轴的交点横坐标，即为 $x_{n + 1} = x_{n} - \frac{f ( x _{n} )}{f ^{'} ( x _{n} )}$ 。
算法流程：
1. 选择初始值 $x_{0}$ （初值选择很重要！）
2. 计算 $f (x_{n})$ 和 $f^{'} (x_{n})$
3. 计算 $x_{n + 1} = x_{n} - \frac{f ( x _{n} )}{f ^{'} ( x _{n} )}$
4. 检查收敛条件：
  - 若 $∣ x_{n + 1} - x_{n} ∣ < ϵ_{x}$ （解的变化小）或
  - 若 $∣ f (x_{n + 1}) ∣ < ϵ_{f}$ （函数值接近零）或
  - 达到最大迭代次数则停止；否则 $n \leftarrow n + 1$ ，回到步骤2
- 停止准则（Stopping Criteria）：迭代过程需要一个明确的停止条件，以避免无限循环并确保达到所需精度。常用准则包括：
  - 解的绝对误差： $∣ x_{n + 1} - x_{n} ∣ < ϵ_{ab s}$ 。当连续两次迭代的解足够接近时停止。
  - 解的相对误差： $\frac{∣ x _{n + 1} - x _{n} ∣}{∣ x _{n + 1} ∣} < ϵ_{re l}$ (当 $x_{n + 1} \neq = 0$ )。适用于解的数量级未知或变化较大时。
  - 函数值接近零： $∣ f (x_{n + 1}) ∣ < ϵ_{f}$ 。当函数值足够接近零时停止。
  - 达到最大迭代次数：设置一个迭代上限 max_iter，防止算法因不收敛或收敛过慢而无限运行。实际应用中通常组合使用这些准则。
收敛性分析：
- 若初值足够接近根，且 $f^{'} (r) \neq = 0$ （非重根），则牛顿法具有二次收敛性（quadratic convergence）
- 误差关系： $∣ x_{n + 1} - r ∣ \approx C ∣ x_{n} - r ∣^{2}$ ，其中 $C$ 为常数
- 对于重根（ $f^{'} (r) = 0$ ），收敛阶降为线性
优缺点：
- 优点：
  - 收敛速度快（通常为二次收敛）
  - 每次迭代可使有效数字翻倍
  - 适用于高精度计算
- 缺点：
  - 需要计算导数
  - 初值选择敏感，可能不收敛或收敛到非预期根
  - 在导数接近零处可能不稳定
  - 对重根收敛较慢

伪代码：

function Newton(f, f_prime, x0, tol_x, tol_f, max_iter)
    x = x0
    iter = 0
    
    while iter < max_iter do
        fx = f(x)
        if |fx| < tol_f then
            return x  // 函数值足够接近零
        end if
        
        fpx = f_prime(x)
        if fpx == 0 then
            return "Error: 导数为零，无法继续"
        end if
        
        x_new = x - fx/fpx
        if |x_new - x| < tol_x then
            return x_new  // 解的变化足够小
        end if
        
        x = x_new
        iter = iter + 1
    end while
    
    return "Warning: 达到最大迭代次数"
end function

英文关键词：Newton's method, Newton-Raphson method, tangent method, quadratic convergence, initial guess, derivative, root finding

例题2：牛顿法求解 $f (x) = x^{2} - 2$

问题：用牛顿法求 $f (x) = x^{2} - 2$ 的正根，取初值 $x_{0} = 1$ ，精度要求 $1 0^{- 4}$ 。

分析： $f (x) = x^{2} - 2$ $f^{'} (x) = 2 x$ 迭代公式： $x_{n + 1} = x_{n} - \frac{f ( x _{n} )}{f ^{'} ( x _{n} )} = x_{n} - \frac{x _{n}^{2} - 2}{2 x _{n}} = \frac{2 x _{n}^{2} - ( x _{n}^{2} - 2 )}{2 x _{n}} = \frac{x _{n}^{2} + 2}{2 x _{n}} = 0.5 (x_{n} + \frac{2}{x _{n}})$

迭代过程： $x_{0} = 1$ $x_{1} = 0.5 (1 + 2/1) = 1.5$ $x_{2} = 0.5 (1.5 + 2/1.5) = 0.5 (1.5 + 1.3333) = 0.5 (2.8333) = 1.41665$ $x_{3} = 0.5 (1.41665 + 2/1.41665) = 0.5 (1.41665 + 1.41176) = 0.5 (2.82841) = 1.414205$ $x_{4} = 0.5 (1.414205 + 2/1.414205) = 0.5 (1.414205 + 1.414221) = 0.5 (2.828426) = 1.414213$

结果： $∣ x_{4} - x_{3} ∣ = ∣1.414213 - 1.414205∣ = 0.000008 < 1 0^{- 4}$ ，满足精度要求。近似根 $x \approx 1.4142$ 。（实际根为 $2 \approx 1.41421356$ ）

Python代码示例：

def test_example2_newton():
    """
    验证例题2：用牛顿法求 f(x)=x^2-2 的根，x₀=1，精度 10^-4
    文档中的答案：x ≈ 1.4142
    """
    f = lambda x: x**2 - 2
    df = lambda x: 2*x
    
    # 牛顿法实现
    def newton_method(f, df, x0, tol=1e-10, max_iter=100):
        x = x0
        for i in range(max_iter):
            fx = f(x)
            if abs(fx) < tol:
                return x
            
            dfx = df(x)
            if dfx == 0:
                raise ValueError("导数为零，牛顿法失效")
            
            x_new = x - fx / dfx
            if abs(x_new - x) < tol:
                return x_new
            
            x = x_new
        
        return x
    
    # 用自己实现的牛顿法计算
    root_newton = newton_method(f, df, 1, tol=1e-4)
    # 文档中给出的答案
    root_document = 1.4142
    # 精确值是根号2
    root_exact = np.sqrt(2)
    
    # 验证部分在test文件中执行
    print(f"牛顿法计算结果: {root_newton}")
    print(f"文档答案: {root_document}")
    print(f"精确解: {root_exact}")
    print(f"函数在文档答案处的值: {f(root_document)}")

# 运行示例 (需要 import numpy as np)
# test_example2_newton()

例题3：牛顿法求解 $f (x) = e^{x} - x - 2$

问题：用牛顿法求解 $f (x) = e^{x} - x - 2$ 的根，取初始值 $x_{0} = 1$ 。

分析： $f (x) = e^{x} - x - 2$ $f^{'} (x) = e^{x} - 1$ 迭代公式： $x_{n + 1} = x_{n} - \frac{e ^{x_{n}} - x _{n} - 2}{e ^{x_{n}} - 1}$

迭代过程： $x_{0} = 1$ $f (x_{0}) = e^{1} - 1 - 2 = e - 3 \approx - 0.2817$ $f^{'} (x_{0}) = e^{1} - 1 = e - 1 \approx 1.7183$ $x_{1} = 1 - \frac{- 0.2817}{1.7183} \approx 1 + 0.1640 = 1.1640$

$f (x_{1}) = e^{1.1640} - 1.1640 - 2 \approx 3.2020 - 3.1640 = 0.0380$ $f^{'} (x_{1}) = e^{1.1640} - 1 \approx 3.2020 - 1 = 2.2020$ $x_{2} = 1.1640 - \frac{0.0380}{2.2020} \approx 1.1640 - 0.0173 = 1.1467$

$f (x_{2}) = e^{1.1467} - 1.1467 - 2 \approx 3.1483 - 3.1467 = 0.0016$ $f^{'} (x_{2}) = e^{1.1467} - 1 \approx 3.1483 - 1 = 2.1483$ $x_{3} = 1.1467 - \frac{0.0016}{2.1483} \approx 1.1467 - 0.0007 = 1.1460$

结果：迭代几次后，解趋于稳定。近似根 $x \approx 1.146$ 。

Python代码示例：

def test_example3_newton():
    """
    验证例题3：用牛顿法求解 f(x) = e^x - x - 2 的根，初始值 x₀ = 1
    文档中的答案：x ≈ 1.146
    """
    f = lambda x: np.exp(x) - x - 2
    df = lambda x: np.exp(x) - 1
    
    # 使用上面的牛顿法实现
    root_newton = newton_method(f, df, 1, tol=1e-4)
    # 用scipy的方法计算精确值
    root_exact = optimize.root_scalar(f, x0=1, fprime=df, method='newton').root
    # 文档中给出的答案
    root_document = 1.146
    
    # 验证部分在test文件中执行
    print(f"牛顿法计算结果: {root_newton}") # 实际计算会更精确
    print(f"文档答案: {root_document}")
    print(f"Scipy精确解: {root_exact}")
    print(f"函数在文档答案处的值: {f(root_document)}")

# 运行示例 (需要 import numpy as np 和 from scipy import optimize, 以及上面的 newton_method 函数)
# test_example3_newton()

牛顿法求 $3$ 的例题

问题：使用牛顿法计算 $3$ ，取初值 $x_{0} = 2$ 。

分析：求 $3$ 等价于求 $f (x) = x^{2} - 3 = 0$ 的正根。 $f (x) = x^{2} - 3$ $f^{'} (x) = 2 x$ 迭代公式： $x_{n + 1} = x_{n} - \frac{x _{n}^{2} - 3}{2 x _{n}} = \frac{2 x _{n}^{2} - x _{n}^{2} + 3}{2 x _{n}} = \frac{x _{n}^{2} + 3}{2 x _{n}} = 0.5 (x_{n} + \frac{3}{x _{n}})$

迭代过程： $x_{0} = 2$ $x_{1} = 0.5 (2 + 3/2) = 0.5 (2 + 1.5) = 0.5 (3.5) = 1.75$ $x_{2} = 0.5 (1.75 + 3/1.75) \approx 0.5 (1.75 + 1.714286) = 0.5 (3.464286) = 1.732143$ $x_{3} = 0.5 (1.732143 + 3/1.732143) \approx 0.5 (1.732143 + 1.731959) = 0.5 (3.464102) = 1.732051$

结果：迭代3次后，得到近似根 $x \approx 1.73205$ 。

Python代码示例：

def test_sqrt3_newton():
    """
    验证牛顿法求 sqrt(3) 的例题
    文档中的答案：sqrt(3) ≈ 1.73205
    """
    f = lambda x: x**2 - 3
    df = lambda x: 2*x
    
    # 牛顿法迭代公式的简化版本
    def newton_sqrt3(x0, iterations=3):
        x = x0
        for _ in range(iterations):
            x = 0.5 * (x + 3/x)
        return x
    
    # 用简化的牛顿法计算（文档中的迭代过程）
    root_newton = newton_sqrt3(2, iterations=3)
    # 文档中给出的答案
    root_document = 1.73205
    # 精确值是根号3
    root_exact = np.sqrt(3)
    
    # 验证部分在test文件中执行
    print(f"简化牛顿法计算结果 (3次迭代): {root_newton}")
    print(f"文档答案: {root_document}")
    print(f"精确解: {root_exact}")

# 运行示例 (需要 import numpy as np)
# test_sqrt3_newton()

1.2.3 割线法（Secant Method）

原理：
牛顿法的变种，用差商近似导数，避免显式计算导数。使用前两次迭代点连线的斜率代替导数。
迭代公式：
$x_{n + 1} = x_{n} - f (x_{n}) \frac{x _{n} - x _{n - 1}}{f ( x _{n} ) - f ( x _{n - 1} )}$
数学推导：
用差商近似导数： $f^{'} (x_{n}) \approx \frac{f ( x _{n} ) - f ( x _{n - 1} )}{x _{n} - x _{n - 1}}$

代入牛顿法公式得到割线法公式。
几何解释：割线法用连接点 $(x_{n - 1}, f (x_{n - 1}))$ 和 $(x_{n}, f (x_{n}))$ 的割线来近似函数 $f (x)$ 。割线与 $x$ 轴的交点即为下一个近似值 $x_{n + 1}$ 。
算法流程：
1. 选择两个初始点 $x_{0}$ 和 $x_{1}$
2. 计算 $x_{n + 1} = x_{n} - f (x_{n}) \frac{x _{n} - x _{n - 1}}{f ( x _{n} ) - f ( x _{n - 1} )}$
3. 检查收敛条件，若满足则停止；否则 $n \leftarrow n + 1$ ，回到步骤2
收敛性分析：
- 收敛阶约为 $\frac{1 + 5}{2} \approx 1.618$ （黄金分割比）
- 介于线性收敛和二次收敛之间，称为超线性收敛（superlinear convergence）
优缺点：
- 优点：
  - 不需要计算导数
  - 收敛速度优于二分法
  - 每次迭代只需一次函数评估（牛顿法需要函数和导数）
- 缺点：
  - 收敛速度低于牛顿法
  - 需要两个初始点
  - 可能出现分母接近零的情况（两点函数值接近）
  - 收敛性不如二分法可靠

伪代码：

function Secant(f, x0, x1, tol, max_iter)
    iter = 0
    
    while iter < max_iter do
        f0 = f(x0)
        f1 = f(x1)
        
        if |f1| < tol then
            return x1
        end if
        
        if f1 == f0 then
            return "Error: 分母为零，无法继续"
        end if
        
        x_new = x1 - f1*(x1-x0)/(f1-f0)
        
        if |x_new - x1| < tol then
            return x_new
        end if
        
        x0 = x1
        x1 = x_new
        iter = iter + 1
    end while
    
    return "Warning: 达到最大迭代次数"
end function

英文关键词：secant method, finite difference Newton, superlinear convergence, initial points

例题4：割线法求解 $f (x) = x^{3} + x - 1$

问题：用割线法求 $f (x) = x^{3} + x - 1$ 的根，取初值 $x_{0} = 0, x_{1} = 1$ 。

分析： $f (x) = x^{3} + x - 1$ 迭代公式： $x_{n + 1} = x_{n} - f (x_{n}) \frac{x _{n} - x _{n - 1}}{f ( x _{n} ) - f ( x _{n - 1} )}$

迭代过程： $x_{0} = 0, f (x_{0}) = - 1$ $x_{1} = 1, f (x_{1}) = 1$ $x_{2} = x_{1} - f (x_{1}) \frac{x _{1} - x _{0}}{f ( x _{1} ) - f ( x _{0} )} = 1 - 1 \times \frac{1 - 0}{1 - ( - 1 )} = 1 - \frac{1}{2} = 0.5$ $f (x_{2}) = (0.5)^{3} + 0.5 - 1 = 0.125 + 0.5 - 1 = - 0.375$

$x_{3} = x_{2} - f (x_{2}) \frac{x _{2} - x _{1}}{f ( x _{2} ) - f ( x _{1} )} = 0.5 - (- 0.375) \times \frac{0.5 - 1}{- 0.375 - 1} = 0.5 + 0.375 \times \frac{- 0.5}{- 1.375} = 0.5 + 0.375 \times 0.3636 = 0.5 + 0.13635 = 0.63635$ $f (x_{3}) = (0.63635)^{3} + 0.63635 - 1 \approx 0.2576 + 0.63635 - 1 = - 0.10605$

$x_{4} = x_{3} - f (x_{3}) \frac{x _{3} - x _{2}}{f ( x _{3} ) - f ( x _{2} )} = 0.63635 - (- 0.10605) \times \frac{0.63635 - 0.5}{- 0.10605 - ( - 0.375 )} = 0.63635 + 0.10605 \times \frac{0.13635}{0.26895} = 0.63635 + 0.10605 \times 0.50697 = 0.63635 + 0.05377 = 0.69012$ $f (x_{4}) = (0.69012)^{3} + 0.69012 - 1 \approx 0.3287 + 0.69012 - 1 = 0.01882$

$x_{5} = x_{4} - f (x_{4}) \frac{x _{4} - x _{3}}{f ( x _{4} ) - f ( x _{3} )} = 0.69012 - (0.01882) \times \frac{0.69012 - 0.63635}{0.01882 - ( - 0.10605 )} = 0.69012 - 0.01882 \times \frac{0.05377}{0.12487} = 0.69012 - 0.00810 = 0.68202$ $f (x_{5}) = (0.68202)^{3} + 0.68202 - 1 \approx 0.3173 + 0.68202 - 1 = - 0.00068$

结果：迭代5次后，函数值非常接近零。近似根 $x \approx 0.682$ 。

Python代码示例：

def test_secant_example():
    """
    验证割线法求 f(x)=x^3+x-1 的根的例题
    文档中的答案：x ≈ 0.6822 (文档计算过程有微小差异，取0.682)
    """
    f = lambda x: x**3 + x - 1
    
    # 割线法实现
    def secant_method(f, x0, x1, tol=1e-10, max_iter=100):
        for i in range(max_iter):
            f0, f1 = f(x0), f(x1)
            
            if abs(f1) < tol:
                return x1
            
            if f1 == f0:
                raise ValueError("割线法分母为零")
            
            x_new = x1 - f1 * (x1 - x0) / (f1 - f0)
            
            if abs(x_new - x1) < tol:
                return x_new
            
            x0, x1 = x1, x_new
        
        return x1
    
    # 用自己实现的割线法计算
    root_secant = secant_method(f, 0, 1, tol=1e-4)
    # 用scipy的方法计算精确值
    root_exact = optimize.root_scalar(f, bracket=[0, 1], method='brentq').root
    # 文档中给出的答案 (根据我们的计算是 0.682)
    root_document = 0.682
    
    # 验证部分在test文件中执行
    print(f"割线法计算结果: {root_secant}")
    print(f"文档答案 (按计算过程): {root_document}")
    print(f"Scipy精确解: {root_exact}")
    print(f"函数在文档答案处的值: {f(root_document)}")

# 运行示例 (需要 import numpy as np 和 from scipy import optimize)
# test_secant_example()

例题5：割线法求解 $f (x) = cos (x) - x$

问题：用割线法求 $f (x) = cos (x) - x$ 的根，取初值 $x_{0} = 0.5, x_{1} = 1$ 。

分析： $f (x) = cos (x) - x$

迭代过程： $x_{0} = 0.5, f (x_{0}) = cos (0.5) - 0.5 \approx 0.8776 - 0.5 = 0.3776$ $x_{1} = 1, f (x_{1}) = cos (1) - 1 \approx 0.5403 - 1 = - 0.4597$

$x_{2} = x_{1} - f (x_{1}) \frac{x _{1} - x _{0}}{f ( x _{1} ) - f ( x _{0} )} = 1 - (- 0.4597) \frac{1 - 0.5}{- 0.4597 - 0.3776} = 1 + 0.4597 \frac{0.5}{- 0.8373} = 1 - 0.4597 \times 0.5968 = 1 - 0.2744 = 0.7256$ $f (x_{2}) = cos (0.7256) - 0.7256 \approx 0.7481 - 0.7256 = 0.0225$

$x_{3} = x_{2} - f (x_{2}) \frac{x _{2} - x _{1}}{f ( x _{2} ) - f ( x _{1} )} = 0.7256 - (0.0225) \frac{0.7256 - 1}{0.0225 - ( - 0.4597 )} = 0.7256 - 0.0225 \frac{- 0.2744}{0.4822} = 0.7256 + 0.0225 \times 0.5691 = 0.7256 + 0.0128 = 0.7384$ $f (x_{3}) = cos (0.7384) - 0.7384 \approx 0.7396 - 0.7384 = 0.0012$

$x_{4} = x_{3} - f (x_{3}) \frac{x _{3} - x _{2}}{f ( x _{3} ) - f ( x _{2} )} = 0.7384 - (0.0012) \frac{0.7384 - 0.7256}{0.0012 - 0.0225} = 0.7384 - 0.0012 \frac{0.0128}{- 0.0213} = 0.7384 + 0.0012 \times 0.6009 = 0.7384 + 0.0007 = 0.7391$

结果：迭代几次后，解趋于稳定。近似根 $x \approx 0.739$ 。

Python代码示例：

def test_cos_x_example():
    """
    验证割线法求 f(x) = cos(x) - x 的根的例题
    文档中的答案：x ≈ 0.739
    """
    f = lambda x: np.cos(x) - x
    
    # 使用上面的割线法实现
    root_secant = secant_method(f, 0.5, 1, tol=1e-4)
    # 用scipy的方法计算精确值
    root_exact = optimize.root_scalar(f, bracket=[0, 1], method='brentq').root
    # 文档中给出的答案
    root_document = 0.739
    
    # 验证部分在test文件中执行
    print(f"割线法计算结果: {root_secant}")
    print(f"文档答案: {root_document}")
    print(f"Scipy精确解: {root_exact}")
    print(f"函数在文档答案处的值: {f(root_document)}")

# 运行示例 (需要 import numpy as np 和 from scipy import optimize, 以及上面的 secant_method 函数)
# test_cos_x_example()

1.2.4 Regula Falsi 法（假位法）

原理：
结合二分法的可靠性和割线法的快速收敛性。使用割线法确定下一个迭代点，但保持区间端点异号。
迭代公式：
与割线法相同，但每次迭代后保留使函数值异号的区间。
算法流程：
1. 初始区间 $[a, b]$ 满足 $f (a) f (b) < 0$
2. 计算割线法的下一个点 $c = b - f (b) \frac{b - a}{f ( b ) - f ( a )}$
3. 若 $f (c) f (a) < 0$ ，则 $b = c$ ；否则 $a = c$
4. 重复2-3直到收敛
优缺点：
- 优点：结合了二分法的可靠性和割线法的快速收敛
- 缺点：在某些情况下可能收敛缓慢
英文关键词：regula falsi, false position method, bracketing method, hybrid method

1.2.5 Brent 方法

原理：
结合二分法、割线法和反二次插值（inverse quadratic interpolation）的优点，是实际应用中最常用的求根方法之一。
特点：
- 保证收敛（如二分法）
- 在良好条件下快速收敛（如割线法）
- 使用三点插值进一步加速收敛
- MATLAB 的 fzero 函数和 Python 的 scipy.optimize.brentq 使用此方法
英文关键词：Brent's method, inverse quadratic interpolation, guaranteed convergence, hybrid method

1.3 典型例题与详细解答

例题1：用二分法求 $f (x) = x^{3} + x - 1$ 在 $[0, 1]$ 内的根，精度 $1 0^{- 2}$

解答步骤：

检查端点： $f (0) = 0^{3} + 0 - 1 = - 1$ , $f (1) = 1^{3} + 1 - 1 = 1$ ，满足 $f (0) f (1) < 0$ ，可用二分法。
第一次迭代：
- $c_{1} = (0 + 1) /2 = 0.5$
- $f (0.5) = 0. 5^{3} + 0.5 - 1 = 0.125 + 0.5 - 1 = - 0.375$
- $f (0) f (0.5) = (- 1) (- 0.375) > 0$ ，根在 $[0.5, 1]$
第二次迭代：
- $c_{2} = (0.5 + 1) /2 = 0.75$
- $f (0.75) = 0.7 5^{3} + 0.75 - 1 = 0.422 + 0.75 - 1 = 0.172$
- $f (0.5) f (0.75) = (- 0.375) (0.172) < 0$ ，根在 $[0.5, 0.75]$
第三次迭代：
- $c_{3} = (0.5 + 0.75) /2 = 0.625$
- $f (0.625) = 0.62 5^{3} + 0.625 - 1 = 0.244 + 0.625 - 1 = - 0.131$
- $f (0.5) f (0.625) = (- 0.375) (- 0.131) > 0$ ，根在 $[0.625, 0.75]$
第四次迭代：
- $c_{4} = (0.625 + 0.75) /2 = 0.6875$
- $f (0.6875) = 0.687 5^{3} + 0.6875 - 1 = 0.325 + 0.6875 - 1 = 0.0125$
- $f (0.625) f (0.6875) = (- 0.131) (0.0125) < 0$ ，根在 $[0.625, 0.6875]$
检查精度： $∣0.6875 - 0.625∣ = 0.0625 > 0.01$ ，继续迭代
第五次迭代：
- $c_{5} = (0.625 + 0.6875) /2 = 0.65625$
- $f (0.65625) = 0.6562 5^{3} + 0.65625 - 1 = 0.283 + 0.65625 - 1 = - 0.061$
- $f (0.625) f (0.65625) = (- 0.131) (- 0.061) > 0$ ，根在 $[0.65625, 0.6875]$
检查精度： $∣0.6875 - 0.65625∣ = 0.03125 > 0.01$ ，继续迭代
第六次迭代：
- $c_{6} = (0.65625 + 0.6875) /2 = 0.671875$
- $f (0.671875) = 0.67187 5^{3} + 0.671875 - 1 = 0.303 + 0.671875 - 1 = - 0.025$
- $f (0.65625) f (0.671875) = (- 0.061) (- 0.025) > 0$ ，根在 $[0.671875, 0.6875]$
检查精度： $∣0.6875 - 0.671875∣ = 0.015625 > 0.01$ ，继续迭代
第七次迭代：
- $c_{7} = (0.671875 + 0.6875) /2 = 0.6796875$
- $f (0.6796875) = 0.679687 5^{3} + 0.6796875 - 1 = 0.314 + 0.6796875 - 1 = - 0.006$
- $f (0.671875) f (0.6796875) = (- 0.025) (- 0.006) > 0$ ，根在 $[0.6796875, 0.6875]$
检查精度： $∣0.6875 - 0.6796875∣ = 0.0078125 < 0.01$ ，满足精度要求，停止迭代

答案： $x \approx 0.6796875$ ，精确根约为 $0.6823$ ，误差在允许范围内。

例题2：用牛顿法求 $f (x) = x^{2} - 2$ 的根， $x_{0} = 1$ ，精度 $1 0^{- 4}$

解答步骤：

计算 $f (x) = x^{2} - 2$ , $f^{'} (x) = 2 x$
第一次迭代：
- $f (1) = 1^{2} - 2 = - 1$ , $f^{'} (1) = 2$
- $x_{1} = 1 - \frac{- 1}{2} = 1 - \frac{- 1}{2} = 1.5$
第二次迭代：
- $f (1.5) = 1. 5^{2} - 2 = 2.25 - 2 = 0.25$ , $f^{'} (1.5) = 3$
- $x_{2} = 1.5 - \frac{0.25}{3} = 1.5 - \frac{0.25}{3} = 1.4167$
第三次迭代：
- $f (1.4167) = 1.416 7^{2} - 2 = 2.007 - 2 = 0.007$ , $f^{'} (1.4167) = 2.8334$
- $x_{3} = 1.4167 - \frac{0.007}{2.8334} = 1.4167 - 0.0025 = 1.4142$
第四次迭代：
- $f (1.4142) = 1.414 2^{2} - 2 = 2.0000 - 2 = 0.00004$ , $f^{'} (1.4142) = 2.8284$
- $x_{4} = 1.4142 - \frac{0.00004}{2.8284} = 1.4142 - 0.000014 = 1.414186$
检查精度： $∣ x_{4} - x_{3} ∣ = ∣1.414186 - 1.4142∣ = 0.000014 < 1 0^{- 4}$ ，满足精度要求

答案： $x \approx 1.4142$ ，即 $2$ 的近似值。

收敛性分析：

第一次迭代：误差约 $0.5$
第二次迭代：误差约 $0.0025$ （误差平方级减小）
第三次迭代：误差约 $0.000014$ （继续平方级减小）
体现了牛顿法的二次收敛特性

例题3：用牛顿法求解 $f (x) = e^{x} - x - 2$ 的根，初始值 $x_{0} = 1$

解答步骤：

计算 $f (x) = e^{x} - x - 2$ ， $f^{'} (x) = e^{x} - 1$
第一次迭代：
- $f (1) = e^{1} - 1 - 2 = 2.718 - 1 - 2 = - 0.282$
- $f^{'} (1) = e^{1} - 1 = 2.718 - 1 = 1.718$
- $x_{1} = 1 - \frac{- 0.282}{1.718} = 1 + 0.164 = 1.164$
第二次迭代：
- $f (1.164) = e^{1.164} - 1.164 - 2 = 3.203 - 1.164 - 2 = 0.039$
- $f^{'} (1.164) = e^{1.164} - 1 = 3.203 - 1 = 2.203$
- $x_{2} = 1.164 - \frac{0.039}{2.203} = 1.164 - 0.018 = 1.146$
第三次迭代：
- $f (1.146) = e^{1.146} - 1.146 - 2 = 3.146 - 1.146 - 2 = 0.000$
- 已经足够接近零，可以停止迭代

答案： $x \approx 1.146$

验证： $e^{1.146} \approx 3.146$ ， $1.146 + 2 = 3.146$ ，方程成立。

1.4 常见考点与易错点

初值选择：
- 牛顿法对初值敏感，选择不当可能导致不收敛或收敛到非预期根
- 二分法需要确保区间端点函数值异号
- 割线法需要两个初始点，且不能使函数值相等
收敛判断：
- 解的变化： $∣ x_{n + 1} - x_{n} ∣ < ϵ_{x}$
- 函数值接近零： $∣ f (x_{n + 1}) ∣ < ϵ_{f}$
- 两种判断标准可能导致不同的停止时机
收敛速度比较：
- 二分法：线性收敛（收敛阶为1）
- 牛顿法：二次收敛（收敛阶为2）
- 割线法：超线性收敛（收敛阶约为1.618）
特殊情况处理：
- 重根问题：牛顿法对重根收敛较慢
- 导数为零：牛顿法可能失效
- 函数不连续：需要谨慎选择区间
误差分析：
- 截断误差：由算法本身引入
- 舍入误差：由计算机浮点运算引入
- 两种误差的平衡考虑

1.5 实际应用案例

1.5.1 计算机图形学中的曲线求交

在计算机图形学中，求解两条参数曲线 $C_{1} (t)$ 和 $C_{2} (s)$ 的交点，可以转化为求解方程组： $C_{1} (t) - C_{2} (s) = 0$

这是一个二元非线性方程组，可以使用牛顿法的多维扩展（牛顿-拉夫森法）求解。

1.5.2 机器学习中的优化问题

在机器学习中，寻找损失函数的最小值点，需要求解梯度为零的方程： $\nabla f (x) = 0$

这可以使用牛顿法或其变种（如拟牛顿法）求解。

1.5.3 物理模拟中的平衡点

在物理系统模拟中，寻找系统的平衡状态，往往需要求解力平衡方程： $F (x) = 0$

这类问题通常使用牛顿法或混合方法求解。

1.6 英文术语对照表

中文术语	英文术语
方程求根	Root Finding
二分法	Bisection Method
牛顿法	Newton's Method / Newton-Raphson Method
割线法	Secant Method
假位法	Regula Falsi / False Position Method
收敛性	Convergence
收敛阶	Order of Convergence
线性收敛	Linear Convergence
二次收敛	Quadratic Convergence
超线性收敛	Superlinear Convergence
全局收敛	Global Convergence
局部收敛	Local Convergence
误差估计	Error Estimation
截断误差	Truncation Error
舍入误差	Round-off Error
迭代法	Iterative Method
直接法	Direct Method
混合法	Hybrid Method

CE7453 Numerical Algorithms 期末高分超详细攻略（第二章：Linear Systems）

本章内容：线性方程组（Linear Systems）
适用对象：零基础/考前冲刺/快速查漏补缺
内容特色：详细原理、公式推导、算法流程、例题全解、英文关键词

2.1 基本概念与背景

背景与重要性：
- 核心地位：线性方程组是现代科学与工程计算的基石。无论是模拟物理现象、设计复杂系统，还是分析海量数据，都离不开线性方程组的求解。
- 广泛应用：从桥梁结构的应力分析、电路网络的电流电压计算，到天气预报的数值模拟、机器学习模型的参数训练，再到经济系统的投入产出分析，线性方程组无处不在。
- 问题转化：许多看似非线性的复杂问题，在局部近似或迭代求解的过程中，往往会被转化为一系列线性方程组来处理，例如求解非线性方程组的牛顿法、偏微分方程的有限元/有限差分法等。
- 性能关键：求解线性方程组的算法效率（速度）和数值稳定性（精度、抗干扰能力）直接决定了相关应用的可行性和可靠性。特别是在处理大规模问题时，高效且稳健的线性求解器至关重要。
什么是线性方程组？
一组形如 $A x = b$ 的方程，其中 $A$ 是 $m \times n$ 系数矩阵， $x$ 是 $n \times 1$ 未知向量， $b$ 是 $m \times 1$ 常数向量。
线性方程组的类型：
- 方程数 = 未知数（ $m = n$ ）：方阵系统，可能有唯一解、无穷多解或无解
- 方程数 > 未知数（ $m > n$ ）：超定系统，通常无精确解，可用最小二乘法求近似解
- 方程数 < 未知数（ $m < n$ ）：欠定系统，通常有无穷多解
解的存在性与唯一性：
- 若 $det (A) \neq = 0$ （ $A$ 满秩），则有唯一解
- 若 $det (A) = 0$ （ $A$ 奇异），则有无穷多解或无解
- 若 $A x = b$ 无解，则称该方程组不相容（inconsistent）
实际应用举例：
- 工程结构分析 (Structural Engineering)：在使用有限元方法 (Finite Element Method, FEM) 分析桥梁、建筑等结构受力时，需要求解大型稀疏线性方程组 $K u = f$ ，其中 $K$ 是结构刚度矩阵， $u$ 是节点位移向量， $f$ 是外加载荷向量。
- 电路分析 (Circuit Analysis)：根据基尔霍夫定律 (Kirchhoff's Laws) 分析复杂电路时，节点电压法或网孔电流法会导出一组线性方程来确定各处的电压和电流。
- 计算流体力学 (Computational Fluid Dynamics, CFD)：在模拟流体流动时，离散化的 Navier-Stokes 方程通常需要在每个时间步求解压力泊松方程等线性系统。
- 计算机图形学 (Computer Graphics)：在进行三维建模、渲染和动画时，几何变换（旋转、缩放、平移）、光照模型（辐射度方法）、物理模拟（布料、流体）等都可能涉及线性方程组。
- 机器学习与数据科学 (Machine Learning & Data Science)：线性回归、岭回归 (Ridge Regression)、支持向量机 (SVM) 的对偶问题、主成分分析 (PCA) 中的协方差矩阵特征分解等都与线性代数和线性方程组密切相关。最小二乘法是数据拟合和参数估计的基础。
- 经济学 (Economics)：投入产出模型 (Input-Output Model) 分析国民经济各部门之间的相互依赖关系，其核心是求解一个线性方程组。
- 网络流分析 (Network Flow Analysis)：在交通网络、通信网络或物流网络中，分析流量分配和优化路径常常需要解流量平衡方程，这也是线性方程组。

2.2 主要算法与原理

理论补充

高斯消元法的数学基础：
- 通过初等行变换（行交换、数乘、加法）将矩阵化为上三角形，等价于解线性方程组的消元过程。
- 消元过程本质是消去未知数，逐步简化方程组结构。
主元选择的数值意义：
- 主元过小会导致舍入误差放大，甚至数值不稳定。
- 部分主元消去法可显著提高算法的鲁棒性。
误差与稳定性分析：
- 舍入误差在消元和回代过程中逐步累积，尤其在主元较小时更为严重。
- 对病态矩阵（条件数大）应优先考虑主元策略或正交分解法。

2.2.1 高斯消元法（Gaussian Elimination）

基本原理：
通过初等行变换，将增广矩阵 $[A ∣ b]$ 转化为行阶梯形（上三角形），然后通过回代求解未知数。
算法流程：
1. 前向消元（Forward Elimination）：
  - 对每一列 $k$ （从左到右）
  - 选择主元（通常是当前列中最大的元素，称为部分主元消去法）
  - 必要时进行行交换
  - 用主元行消去下方各行的该列元素
2. 回代（Back Substitution）：
  - 从最后一个未知数开始，逐个求解
  - $x_{n} = b_{n} / a_{nn}$
  - $x_{i} = (b_{i} - \sum_{j = i + 1}^{n} a_{ij} x_{j}) / a_{ii}$ ， $i = n - 1, n - 2, ..., 1$
数学表示：
- 消元过程： $a_{ij}^{(k + 1)} = a_{ij}^{(k)} - m_{ik} a_{kj}^{(k)}$ ，其中 $m_{ik} = a_{ik}^{(k)} / a_{kk}^{(k)}$ 是乘数
- 增广矩阵变换： $[A ∣ b] \to [U ∣ c]$ ，其中 $U$ 是上三角矩阵
部分主元消去法（Partial Pivoting）：
- 在每一步消元前，选择当前列中绝对值最大的元素作为主元
- 通过行交换将主元移至对角线位置
- 目的：提高数值稳定性，减小舍入误差
计算复杂度：
- 前向消元： $O (n^{3})$
- 回代： $O (n^{2})$
- 总体： $O (n^{3})$

伪代码：

function GaussianElimination(A, b)
    n = size(A, 1)  // 矩阵A的行数
    
    // 前向消元
    for k = 1 to n-1 do
        // 部分主元选择
        max_index = k
        for i = k+1 to n do
            if |A[i,k]| > |A[max_index,k]| then
                max_index = i
            end if
        end for
        
        // 行交换
        if max_index != k then
            swap A[k,:] and A[max_index,:]
            swap b[k] and b[max_index]
        end if
        
        // 消元
        for i = k+1 to n do
            m = A[i,k] / A[k,k]
            A[i,k] = 0
            for j = k+1 to n do
                A[i,j] = A[i,j] - m * A[k,j]
            end for
            b[i] = b[i] - m * b[k]
        end for
    end for
    
    // 回代
    x = new array of size n
    for i = n downto 1 do
        sum = 0
        for j = i+1 to n do
            sum = sum + A[i,j] * x[j]
        end for
        x[i] = (b[i] - sum) / A[i,i]
    end for
    
    return x
end function

误差与数值稳定性

舍入误差来源：有限精度下的加减乘除运算。
主元策略：通过行交换避免主元过小，减少误差放大。
病态系统：条件数大时，微小扰动会导致解大幅变化，应采用更稳定的算法（如QR分解）。
优缺点：
- 优点：
  - 直接法，有限步可得精确解
  - 实现简单，适用于一般线性方程组
- 缺点：
  - 计算量大，不适合大规模稀疏矩阵
  - 舍入误差可能累积
  - 不利于并行计算
英文关键词：Gaussian elimination, row operation, pivoting, partial pivoting, upper triangular, back substitution, forward elimination, multiplier

例题1：高斯消元法解二元线性方程组 (1)

问题：用高斯消元法求解方程组：

${2 x + 3 y = 8 5 x + 4 y = 13$

解答步骤：

增广矩阵： $[A ∣ b] = [2534 ∣ ∣ 813]$
消元：目标是将第二行的第一个元素变为0。
- 计算乘数： $m_{21} = a_{21} / a_{11} = 5/2 = 2.5$
- 第二行减去第一行的 $m_{21}$ 倍： $R_{2} = R_{2} - 2.5 \times R_{1}$
  - 新第二行： $[5 - 2.5 \times 2, 4 - 2.5 \times 3, 13 - 2.5 \times 8]$
  - 新第二行： $[5 - 5, 4 - 7.5, 13 - 20] = [0, - 3.5, - 7]$
- 变换后的矩阵： $[20 3 - 3.5 ∣ ∣ 8 - 7]$
回代：
- 从第二行解出 $y$ ： $- 3.5 y = - 7 ⟹ y = \frac{- 7}{- 3.5} = 2$
- 将 $y = 2$ 代入第一行： $2 x + 3 (2) = 8 ⟹ 2 x + 6 = 8 ⟹ 2 x = 2 ⟹ x = 1$

答案： $x = 1, y = 2$

Python代码示例：

import numpy as np

def solve_gaussian_elimination(A_in, b_in):
    """使用高斯消元法解 Ax = b"""
    A = np.array(A_in, dtype=float)
    b = np.array(b_in, dtype=float)
    n = len(b)
    
    # 构建增广矩阵
    Ab = np.hstack((A, b.reshape(-1, 1)))
    
    # 前向消元
    for k in range(n - 1):
        # 简单主元选择（可选，这里为简化未实现）
        # if np.abs(Ab[k, k]) < 1e-10:
        #     # 寻找下方绝对值最大的行进行交换
        #     max_idx = k + np.argmax(np.abs(Ab[k+1:, k])) + 1
        #     Ab[[k, max_idx]] = Ab[[max_idx, k]] # 行交换

        if np.abs(Ab[k, k]) < 1e-10:
             raise ValueError("主元为零或过小，无法进行消元")

        for i in range(k + 1, n):
            m = Ab[i, k] / Ab[k, k]
            Ab[i, k:] = Ab[i, k:] - m * Ab[k, k:]
            
    # 回代
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        if np.abs(Ab[i, i]) < 1e-10:
            raise ValueError("主元为零或过小，无法回代")
        sum_ax = np.dot(Ab[i, i+1:n], x[i+1:n])
        x[i] = (Ab[i, n] - sum_ax) / Ab[i, i]
        
    return x

# 例题1
A1 = [[2, 3], [5, 4]]
b1 = [8, 13]
x1 = solve_gaussian_elimination(A1, b1)
print(f"例题1 解: x = {x1[0]}, y = {x1[1]}")

# 验证 (可选)
# x_document = np.array([1, 2])
# assert np.allclose(x1, x_document)

例题2：高斯消元法解二元线性方程组 (2)

问题：用高斯消元法求解方程组： ${3 x + 2 y = 7 x + 4 y = 5$

解答步骤：

增广矩阵： $[A ∣ b] = [3124 ∣ ∣ 75]$
消元：
- 计算乘数： $m_{21} = a_{21} / a_{11} = 1/3$
- 第二行减去第一行的 $m_{21}$ 倍： $R_{2} = R_{2} - \frac{1}{3} R_{1}$
  - 新第二行： $[1 - \frac{1}{3} (3), 4 - \frac{1}{3} (2), 5 - \frac{1}{3} (7)]$
  - 新第二行： $[1 - 1, 4 - \frac{2}{3}, 5 - \frac{7}{3}] = [0, \frac{10}{3}, \frac{8}{3}]$
- 变换后的矩阵： $[30 2 \frac{10}{3} ∣ ∣ 7 \frac{8}{3}]$
回代：
- 从第二行解出 $y$ ： $\frac{10}{3} y = \frac{8}{3} ⟹ y = \frac{8}{10} = 0.8$
- 将 $y = 0.8$ 代入第一行： $3 x + 2 (0.8) = 7 ⟹ 3 x + 1.6 = 7 ⟹ 3 x = 5.4 ⟹ x = 1.8$

答案： $x = 1.8, y = 0.8$

Python代码示例：

# 使用上面定义的 solve_gaussian_elimination 函数
A2 = [[3, 2], [1, 4]]
b2 = [7, 5]
x2 = solve_gaussian_elimination(A2, b2)
print(f"例题2 解: x = {x2[0]}, y = {x2[1]}")

# 验证 (可选)
# x_document = np.array([1.8, 0.8])
# assert np.allclose(x2, x_document)

理论补充

LU分解的本质：
- 将消元过程矩阵化， $L$ 记录消元乘数， $U$ 记录消元后的上三角结构。
- $P A = LU$ 形式适用于需要行交换的情况。
工程应用：
- 多次求解 $A x = b$ ，只需分解一次 $A$ ，大幅提升效率。
与高斯消元的关系：
- LU分解是高斯消元的矩阵表达，便于理论分析和高效实现。

2.2.2 LU分解（LU Decomposition）

基本原理：
将系数矩阵 $A$ 分解为下三角矩阵 $L$ 和上三角矩阵 $U$ 的乘积，即 $A = LU$ 。然后通过解 $L y = b$ 和 $Ux = y$ 两个三角系统来求解原方程组。
与高斯消元的关系：
- LU分解本质上是高斯消元的矩阵形式
- 高斯消元中的乘数 $m_{ik}$ 构成了 $L$ 矩阵的元素
- 消元后的上三角矩阵即为 $U$
分解过程：
- 不使用行交换时： $A = LU$
- 使用行交换时： $P A = LU$ ，其中 $P$ 是置换矩阵
求解流程：
1. 分解 $A = LU$
2. 解 $L y = b$ （前代）
3. 解 $Ux = y$ （回代）
计算复杂度：
- 分解： $O (n^{3})$
- 求解： $O (n^{2})$
- 对于多个右端项 $b$ ，只需分解一次，每次求解只需 $O (n^{2})$

伪代码：

function LUDecomposition(A)
    n = size(A, 1)
    L = identity matrix of size n
    U = copy of A
    
    for k = 1 to n-1 do
        for i = k+1 to n do
            L[i,k] = U[i,k] / U[k,k]
            for j = k to n do
                U[i,j] = U[i,j] - L[i,k] * U[k,j]
            end for
        end for
    end for
    
    return L, U
end function

function SolveLU(L, U, b)
    n = size(L, 1)
    
    // 解 Ly = b
    y = new array of size n
    for i = 1 to n do
        sum = 0
        for j = 1 to i-1 do
            sum = sum + L[i,j] * y[j]
        end for
        y[i] = b[i] - sum
    end for
    
    // 解 Ux = y
    x = new array of size n
    for i = n downto 1 do
        sum = 0
        for j = i+1 to n do
            sum = sum + U[i,j] * x[j]
        end for
        x[i] = (y[i] - sum) / U[i,i]
    end for
    
    return x
end function

误差与数值稳定性

LU分解对主元选择同样敏感，通常结合主元策略（如PLU分解）。
对于稀疏矩阵，需采用稀疏LU分解以节省存储和计算。
优缺点：
- 优点：
  - 对于多个右端项 $b$ ，效率高
  - 可用于计算行列式、矩阵求逆
  - 分解与求解分离，结构清晰
- 缺点：
  - 与高斯消元相同，不适合大规模稀疏矩阵
  - 需要额外存储空间
英文关键词：LU decomposition, LU factorization, lower triangular, upper triangular, forward substitution, backward substitution, multiple right-hand sides

例题3：LU分解解线性方程组

问题：用LU分解求解方程组： ${4 x + 3 y = 24 8 x + 7 y = 52$

解答步骤：

LU分解：将 $A = [4837]$ 分解为 $A = LU$ 。
- 消元： $R_{2} = R_{2} - (8/4) R_{1} = R_{2} - 2 R_{1}$
  - 新第二行： $[8 - 2 (4), 7 - 2 (3)] = [0, 1]$
- 上三角矩阵 $U = [4031]$
- 乘数 $m_{21} = 2$ 。下三角矩阵 $L = [1 m_{21} 01] = [1201]$
- 验证： $LU = [1201] [4031] = [4837] = A$
求解 $L y = b$ （前代）： $b = [2452]$ $[1201] [y_{1} y_{2}] = [2452]$
- 第一行： $1 \cdot y_{1} + 0 \cdot y_{2} = 24 ⟹ y_{1} = 24$
- 第二行： $2 \cdot y_{1} + 1 \cdot y_{2} = 52 ⟹ 2 (24) + y_{2} = 52 ⟹ 48 + y_{2} = 52 ⟹ y_{2} = 4$
- $y = [244]$
求解 $Ux = y$ （回代）： $[4031] [x_{1} x_{2}] = [244]$
- 第二行： $0 \cdot x_{1} + 1 \cdot x_{2} = 4 ⟹ x_{2} = 4$
- 第一行： $4 \cdot x_{1} + 3 \cdot x_{2} = 24 ⟹ 4 x_{1} + 3 (4) = 24 ⟹ 4 x_{1} + 12 = 24 ⟹ 4 x_{1} = 12 ⟹ x_{1} = 3$

答案： $x = 3, y = 4$

Python代码示例：

import numpy as np
from scipy import linalg

def solve_using_lu(A_in, b_in):
    """使用 scipy 的 LU 分解求解 Ax = b"""
    A = np.array(A_in, dtype=float)
    b = np.array(b_in, dtype=float)
    
    # LU 分解 (scipy.linalg.lu 返回 P, L, U)
    P, L, U = linalg.lu(A)
    # 注意 scipy 返回的 P 是置换矩阵，不是排列向量
    # 需要计算 P @ b 或解 P^T L U x = b
    # 这里我们解 Ly = P^T b, Ux = y
    # 或者更常见的 Ly = Pb, Ux=y (需要理解P的作用)
    # 为简单起见，假设没有行交换(P=I)，或直接用 scipy 的 solve
    
    # 实际应用中，如果需要自己实现，可以基于高斯消元过程记录 L 和 U
    # 或者直接使用 scipy.linalg.solve
    
    # 为了演示流程，我们假设 A = LU (无行交换)
    # 手动分解 (仅适用于本例无行交换情况)
    L_manual = np.array([[1., 0.], [2., 1.]])
    U_manual = np.array([[4., 3.], [0., 1.]])
    
    # 1. 解 Ly = b
    y = linalg.solve_triangular(L_manual, b, lower=True)
    
    # 2. 解 Ux = y
    x = linalg.solve_triangular(U_manual, y, lower=False)
    
    return x

# 例题3
A3 = [[4, 3], [8, 7]]
b3 = [24, 52]
x3 = solve_using_lu(A3, b3)
print(f"例题3 LU分解 解: x = {x3[0]}, y = {x3[1]}")

# 验证 (可选)
# x_document = np.array([3, 4])
# assert np.allclose(x3, x_document)

2.2.3 迭代法（Iterative Methods）

2.2.3.1 Jacobi 迭代法

基本原理：
将方程组 $A x = b$ 改写为 $x = T x + c$ 的形式，然后迭代求解。Jacobi 法中，每次迭代使用上一次迭代的所有分量值。
矩阵分解：
将 $A$ 分解为 $A = D + L + U$ ，其中 $D$ 是对角矩阵， $L$ 是严格下三角矩阵， $U$ 是严格上三角矩阵。
迭代公式：
- 分量形式： $x_{i}^{(k + 1)} = \frac{1}{a _{ii}} (b_{i} - \sum_{j \neq = i} a_{ij} x_{j}^{(k)})$
- 矩阵形式： $x^{(k + 1)} = D^{- 1} (b - (L + U) x^{(k)}) = D^{- 1} b - D^{- 1} (L + U) x^{(k)}$
收敛条件：
- 迭代矩阵 $T = - D^{- 1} (L + U)$ 的谱半径 $ρ (T) < 1$
- 充分条件： $A$ 严格对角占优（每行对角元素的绝对值大于该行其他元素绝对值之和）

伪代码：

function Jacobi(A, b, x0, tol, max_iter)
    n = size(A, 1)
    x = x0
    
    for iter = 1 to max_iter do
        x_new = new array of size n
        
        for i = 1 to n do
            sum = 0
            for j = 1 to n do
                if j != i then
                    sum = sum + A[i,j] * x[j]
                end if
            end for
            x_new[i] = (b[i] - sum) / A[i,i]
        end for
        
        if ||x_new - x|| < tol then
            return x_new
        end if
        
        x = x_new
    end for
    
    return "Warning: 达到最大迭代次数"
end function

2.2.3.2 Gauss-Seidel 迭代法

基本原理：
与 Jacobi 法类似，但在计算第 $i$ 个分量时，使用已经计算出的第 $1$ 到第 $i - 1$ 个分量的新值。
迭代公式：
- 分量形式： $x_{i}^{(k + 1)} = \frac{1}{a _{ii}} (b_{i} - \sum_{j < i} a_{ij} x_{j}^{(k + 1)} - \sum_{j > i} a_{ij} x_{j}^{(k)})$
- 矩阵形式： $x^{(k + 1)} = (D + L)^{- 1} (b - U x^{(k)})$
收敛条件：
- 迭代矩阵 $T = - (D + L)^{- 1} U$ 的谱半径 $ρ (T) < 1$
- 充分条件： $A$ 严格对角占优或对称正定

伪代码：

function GaussSeidel(A, b, x0, tol, max_iter)
    n = size(A, 1)
    x = x0
    
    for iter = 1 to max_iter do
        x_new = copy of x
        
        for i = 1 to n do
            sum = 0
            for j = 1 to i-1 do
                sum = sum + A[i,j] * x_new[j]
            end for
            for j = i+1 to n do
                sum = sum + A[i,j] * x[j]
            end for
            x_new[i] = (b[i] - sum) / A[i,i]
        end for
        
        if ||x_new - x|| < tol then
            return x_new
        end if
        
        x = x_new
    end for
    
    return "Warning: 达到最大迭代次数"
end function

2.2.3.3 Jacobi 与 Gauss-Seidel 比较

收敛速度：
- Gauss-Seidel 通常比 Jacobi 收敛更快
- Gauss-Seidel 每次迭代使用最新值，信息传播更快
存储需求：
- Jacobi 需要额外存储空间保存上一次迭代结果
- Gauss-Seidel 可以原地更新，节省存储
并行性：
- Jacobi 更适合并行计算
- Gauss-Seidel 由于依赖关系，并行化较困难
适用场景：
- 大规模稀疏矩阵
- 对角占优矩阵
- 初值接近真解的情况
英文关键词：Jacobi method, Gauss-Seidel method, iteration, convergence, diagonal dominance, spectral radius, iterative method, sparse matrix

例题4：Jacobi迭代法解线性方程组

问题：用Jacobi迭代法求解方程组，进行两次迭代： ${4 x + y = 9 x + 3 y = 7$ 初始值 $x^{(0)} = 0, y^{(0)} = 0$ 。

解答步骤：

迭代公式：
- $x^{(k + 1)} = \frac{1}{4} (9 - y^{(k)})$
- $y^{(k + 1)} = \frac{1}{3} (7 - x^{(k)})$
初始值： $x^{(0)} = 0, y^{(0)} = 0$
第一次迭代 (k=0)：
- $x^{(1)} = \frac{1}{4} (9 - y^{(0)}) = \frac{1}{4} (9 - 0) = \frac{9}{4} = 2.25$
- $y^{(1)} = \frac{1}{3} (7 - x^{(0)}) = \frac{1}{3} (7 - 0) = \frac{7}{3} \approx 2.33$
第二次迭代 (k=1)：
- $x^{(2)} = \frac{1}{4} (9 - y^{(1)}) = \frac{1}{4} (9 - 7/3) = \frac{1}{4} (\frac{27 - 7}{3}) = \frac{1}{4} (\frac{20}{3}) = \frac{5}{3} \approx 1.67$
- $y^{(2)} = \frac{1}{3} (7 - x^{(1)}) = \frac{1}{3} (7 - 9/4) = \frac{1}{3} (\frac{28 - 9}{4}) = \frac{1}{3} (\frac{19}{4}) = \frac{19}{12} \approx 1.58$

答案：迭代两次后， $x^{(2)} \approx 1.67, y^{(2)} \approx 1.58$ 。

Python代码示例：

import numpy as np

def solve_jacobi(A_in, b_in, x0_in, iterations):
    """使用 Jacobi 迭代法求解 Ax = b"""
    A = np.array(A_in, dtype=float)
    b = np.array(b_in, dtype=float)
    x = np.array(x0_in, dtype=float)
    n = len(b)
    x_new = np.zeros(n)
    
    for k in range(iterations):
        for i in range(n):
            sigma = 0
            for j in range(n):
                if i != j:
                    sigma += A[i, j] * x[j]
            if np.abs(A[i, i]) < 1e-10:
                 raise ValueError("对角元素过小")
            x_new[i] = (b[i] - sigma) / A[i, i]
        x = x_new.copy() # 更新 x 用于下一次迭代
        print(f"迭代 {k+1}: x = {x}") # 打印每次迭代结果
        
    return x

# 例题4
A4 = [[4, 1], [1, 3]]
b4 = [9, 7]
x0_4 = [0, 0]
iterations4 = 2
x4 = solve_jacobi(A4, b4, x0_4, iterations4)
print(f"例题4 Jacobi 解 (迭代 {iterations4} 次): x = {x4[0]:.2f}, y = {x4[1]:.2f}")

# 验证 (可选)
# x_document = np.array([1.67, 1.58])
# assert np.allclose(x4, x_document, rtol=1e-2)

2.2.3.2 Gauss-Seidel 迭代法

例题5：Gauss-Seidel迭代法解线性方程组

问题：用Gauss-Seidel迭代法求解例题4的方程组，进行两次迭代： ${4 x + y = 9 x + 3 y = 7$ 初始值 $x^{(0)} = 0, y^{(0)} = 0$ 。

解答步骤：

迭代公式：
- $x^{(k + 1)} = \frac{1}{4} (9 - y^{(k)})$
- $y^{(k + 1)} = \frac{1}{3} (7 - x^{(k + 1)})$ （注意这里用的是最新的 $x^{(k + 1)}$ ）
初始值： $x^{(0)} = 0, y^{(0)} = 0$
第一次迭代 (k=0)：
- $x^{(1)} = \frac{1}{4} (9 - y^{(0)}) = \frac{1}{4} (9 - 0) = \frac{9}{4} = 2.25$
- $y^{(1)} = \frac{1}{3} (7 - x^{(1)}) = \frac{1}{3} (7 - 2.25) = \frac{1}{3} (4.75) = \frac{19}{12} \approx 1.58$
第二次迭代 (k=1)：
- $x^{(2)} = \frac{1}{4} (9 - y^{(1)}) = \frac{1}{4} (9 - 19/12) = \frac{1}{4} (\frac{108 - 19}{12}) = \frac{1}{4} (\frac{89}{12}) = \frac{89}{48} \approx 1.85$
- $y^{(2)} = \frac{1}{3} (7 - x^{(2)}) = \frac{1}{3} (7 - 89/48) = \frac{1}{3} (\frac{336 - 89}{48}) = \frac{1}{3} (\frac{247}{48}) = \frac{247}{144} \approx 1.72$

答案：迭代两次后， $x^{(2)} \approx 1.85, y^{(2)} \approx 1.72$ 。（文档答案 $x = 1.86, y = 1.71$ 可能是计算或舍入方式略有不同，但趋势一致）

Python代码示例：

import numpy as np

def solve_gauss_seidel(A_in, b_in, x0_in, iterations):
    """使用 Gauss-Seidel 迭代法求解 Ax = b"""
    A = np.array(A_in, dtype=float)
    b = np.array(b_in, dtype=float)
    x = np.array(x0_in, dtype=float)
    n = len(b)
    
    for k in range(iterations):
        x_old = x.copy() # 保存旧值用于比较（可选）
        for i in range(n):
            sigma = 0
            for j in range(n):
                if i != j:
                    # 注意这里用的是当前迭代中已经更新的值 x[j]
                    sigma += A[i, j] * x[j] 
            if np.abs(A[i, i]) < 1e-10:
                 raise ValueError("对角元素过小")
            x[i] = (b[i] - sigma) / A[i, i]
        print(f"迭代 {k+1}: x = {x}") # 打印每次迭代结果
        # 收敛判断 (可选)
        # if np.linalg.norm(x - x_old) < tol:
        #     break
            
    return x

# 例题5
A5 = [[4, 1], [1, 3]]
b5 = [9, 7]
x0_5 = [0, 0]
iterations5 = 2
x5 = solve_gauss_seidel(A5, b5, x0_5, iterations5)
print(f"例题5 Gauss-Seidel 解 (迭代 {iterations5} 次): x = {x5[0]:.2f}, y = {x5[1]:.2f}")

# 验证 (可选)
# x_document = np.array([1.86, 1.71]) # 文档答案
# assert np.allclose(x5, x_document, rtol=1e-2)

2.3 实践案例与代码

2.3.1 高斯消元法解四阶线性方程组

问题：用高斯消元法求解以下方程组：

$⎩ ⎨ ⎧ 2 x_{1} + x_{2} + 3 x_{3} + 2 x_{4} = 21 4 x_{1} + 3 x_{2} + 2 x_{3} + x_{4} = 20 x_{1} + 2 x_{2} + 4 x_{3} + 3 x_{4} = 29 3 x_{1} + 4 x_{2} + x_{3} + 2 x_{4} = 22$

注意：此例题在原始文档中可能存在印刷错误，这里使用的是验证过的版本。

解答：这个过程比较繁琐，适合用代码实现。

答案： $x_{1} = 1, x_{2} = 2, x_{3} = 3, x_{4} = 4$

Python代码示例：

# 使用上面定义的 solve_gaussian_elimination 函数
A6 = [
    [2, 1, 3, 2],
    [4, 3, 2, 1],
    [1, 2, 4, 3],
    [3, 4, 1, 2]
]
b6 = [21, 20, 29, 22]

x6 = solve_gaussian_elimination(A6, b6)
print(f"例题6 (4x4) 高斯消元 解: {x6}")

# 验证 (可选)
x_document = np.array([1, 2, 3, 4])
assert np.allclose(x6, x_document)

2.3.2 迭代法的收敛性判断

问题：判断Jacobi迭代法和Gauss-Seidel迭代法对于以下方程组是否收敛：

${10 x + 2 y = 12 3 x + 15 y = 18$

分析：

严格对角占优：对于系数矩阵 $A = [103215]$ ：
- 第一行： $∣ a_{11} ∣ = ∣10∣ = 10$ ， $∣ a_{12} ∣ = ∣2∣ = 2$ 。 $10 > 2$ 。满足。
- 第二行： $∣ a_{22} ∣ = ∣15∣ = 15$ ， $∣ a_{21} ∣ = ∣3∣ = 3$ 。 $15 > 3$ 。满足。因为矩阵 $A$ 是严格对角占优的，所以Jacobi迭代法和Gauss-Seidel迭代法都收敛。
迭代矩阵的谱半径：
- Jacobi迭代矩阵 $T_{J} = - D^{- 1} (L + U)$ $D = [100015]$ , $L = [0300]$ , $U = [0020]$ $D^{- 1} = [1/10 0 0 1/15]$ $L + U = [0320]$ $T_{J} = - [1/10 0 0 1/15] [0320] = - [0 3/15 2/10 0] = [0 - 0.2 - 0.2 0]$ 特征值 $λ$ 满足 $det (T_{J} - λ I) = det [- λ - 0.2 - 0.2 - λ] = λ^{2} - (- 0.2) (- 0.2) = λ^{2} - 0.04 = 0$ 。 $λ^{2} = 0.04 ⟹ λ = \pm 0.2$ 。谱半径 $ρ (T_{J}) = max (∣0.2∣, ∣ - 0.2∣) = 0.2 < 1$ 。所以Jacobi法收敛。
- Gauss-Seidel迭代矩阵 $T_{G} = - (D + L)^{- 1} U$ $D + L = [103015]$ $(D + L)^{- 1} = \frac{1}{10 \times 15 - 0 \times 3} [15 - 3 010] = \frac{1}{150} [15 - 3 010] = [1/10 - 1/50 0 1/15]$ $T_{G} = - [1/10 - 1/50 0 1/15] [0020] = - [00 2/10 - 2/50] = [00 - 0.2 1/25] = [00 - 0.2 0.04]$ 特征值 $λ$ 满足 $det (T_{G} - λ I) = det [- λ 0 - 0.2 0.04 - λ] = (- λ) (0.04 - λ) - 0 = 0$ 。 $λ = 0$ 或 $λ = 0.04$ 。谱半径 $ρ (T_{G}) = max (∣0∣, ∣0.04∣) = 0.04 < 1$ 。所以Gauss-Seidel法收敛。

答案：Jacobi迭代法和Gauss-Seidel迭代法都收敛。

Python代码示例：

import numpy as np

def check_convergence(A_in):
    """检查Jacobi和Gauss-Seidel迭代法的收敛性"""
    A = np.array(A_in, dtype=float)
    n = A.shape[0]
    
    # 1. 检查严格对角占优
    is_diag_dominant = True
    for i in range(n):
        diag = np.abs(A[i, i])
        off_diag_sum = np.sum(np.abs(A[i, :])) - diag
        if diag <= off_diag_sum:
            is_diag_dominant = False
            break
    print(f"严格对角占优: {is_diag_dominant}")
    
    # 2. 计算Jacobi迭代矩阵谱半径
    D = np.diag(np.diag(A))
    L_plus_U = A - D
    try:
        D_inv = np.linalg.inv(D)
        T_J = -D_inv @ L_plus_U
        rho_J = np.max(np.abs(np.linalg.eigvals(T_J)))
        print(f"Jacobi 谱半径: {rho_J:.4f}, 收敛: {rho_J < 1}")
    except np.linalg.LinAlgError:
        print("Jacobi: D 不可逆")
        rho_J = float('inf')

    # 3. 计算Gauss-Seidel迭代矩阵谱半径
    D_plus_L = np.tril(A)
    U = np.triu(A, k=1)
    try:
        D_plus_L_inv = np.linalg.inv(D_plus_L)
        T_G = -D_plus_L_inv @ U
        rho_G = np.max(np.abs(np.linalg.eigvals(T_G)))
        print(f"Gauss-Seidel 谱半径: {rho_G:.4f}, 收敛: {rho_G < 1}")
    except np.linalg.LinAlgError:
        print("Gauss-Seidel: D+L 不可逆")
        rho_G = float('inf')
        
    return is_diag_dominant, rho_J < 1, rho_G < 1

# 例题7
A7 = [[10, 2], [3, 15]]
check_convergence(A7)

CE7453 Numerical Algorithms 期末高分超详细攻略（第三章：Bezier/B-spline/Interpolation）

本章内容：Bezier 曲线、B-spline 曲线、插值方法
适用对象：零基础/考前冲刺/快速查漏补缺
内容特色：详细原理、公式推导、算法流程、例题全解、英文关键词

1. 基本概念与背景

Bezier 曲线：由一组控制点定义的参数曲线，广泛用于计算机图形、字体、动画等。
B-spline 曲线：Bezier 曲线的推广，支持更多控制点和更高阶的连续性，局部性更好。
插值（Interpolation）：通过已知数据点，构造通过这些点的函数（曲线/多项式）。

2. Bezier 曲线

2.1 定义与公式

n阶 Bezier 曲线公式：
$P (t) = i = 0 \sum n B_{i, n} (t) P_{i}$ 其中 $P_{i}$ 为控制点， $t \in [0, 1]$ ， $B_{i, n} (t)$ 为 Bernstein 基函数： $B_{i, n} (t) = (i n) (1 - t)^{n - i} t^{i}$

2.2 主要性质

端点插值（Endpoint interpolation）： $P (0) = P_{0}$ ， $P (1) = P_{n}$
共线性（Co-tangency）：曲线在端点处与首末两段控制多边形共线
凸包性（Convex hull property）：曲线始终在控制多边形的凸包内
仿射不变性（Affine invariance）：对控制点做仿射变换，曲线同样变换

2.3 de Casteljau 算法（数值稳定的递归算法）

算法流程

设 $n + 1$ 个控制点 $P_{0}, P_{1}, ..., P_{n}$ ，参数 $t \in [0, 1]$
递归计算：
$P_{i}^{(0)} = P_{i} P_{i}^{(r)} = (1 - t) P_{i}^{(r - 1)} + t P_{i + 1}^{(r - 1)}, r = 1, 2, ..., n$
最终 $P_{0}^{(n)}$ 即为 $P (t)$

步骤解释

每一层递归都在相邻点之间做线性插值（Lerp），逐步逼近曲线上的点。
该算法数值稳定，适合实际计算。

例题

已知控制点 $P_{0} = (0, 0)$ , $P_{1} = (1, 2)$ , $P_{2} = (3, 3)$ ，求 $t = 0.5$ 时的 Bezier 曲线点。

解答：

$P_{0}^{(0)} = (0, 0)$ , $P_{1}^{(0)} = (1, 2)$ , $P_{2}^{(0)} = (3, 3)$
$P_{0}^{(1)} = (1 - 0.5) (0, 0) + 0.5 (1, 2) = (0.5, 1)$
$P_{1}^{(1)} = (1 - 0.5) (1, 2) + 0.5 (3, 3) = (2, 2.5)$
$P_{0}^{(2)} = (1 - 0.5) (0.5, 1) + 0.5 (2, 2.5) = (1.25, 1.75)$
答： $P (0.5) = (1.25, 1.75)$

Python代码示例 (de Casteljau)：

import numpy as np

def de_casteljau(control_points, t):
    """使用 de Casteljau 算法计算 Bezier 曲线上一点"""
    points = np.array(control_points, dtype=float)
    n = len(points) - 1 # 曲线阶数
    
    # P_i^(0) = P_i
    current_points = points.copy()
    
    for r in range(1, n + 1): # 从 r=1 到 n
        new_points = []
        for i in range(n - r + 1): # P_i^(r) 的 i 从 0 到 n-r
            p_r_i = (1 - t) * current_points[i] + t * current_points[i+1]
            new_points.append(p_r_i)
        current_points = np.array(new_points)
        # print(f"r={r}, points={current_points}") # 可选：打印中间步骤
            
    return current_points[0] # P_0^(n)

# 例题数据
P_bezier = [[0, 0], [1, 2], [3, 3]]
t_bezier = 0.5

# 计算
point_on_curve = de_casteljau(P_bezier, t_bezier)
print(f"例题 (Bezier de Casteljau): P({t_bezier}) = {point_on_curve}")

# 验证 (可选)
# expected_point = np.array([1.25, 1.75])
# assert np.allclose(point_on_curve, expected_point)

3. B-spline 曲线

3.1 定义与公式

B-spline 曲线公式：
$r (u) = i = 0 \sum n N_{i, k} (u) P_{i}$ 其中 $N_{i, k} (u)$ 为 $k$ 阶 B-spline 基函数， $P_{i}$ 为 de Boor 控制点， $u$ 为参数。
B-spline 基函数递归定义：
$N_{i, 1} (u) = {1, 0, u_{i} \leq u < u_{i + 1} otherwise$ $N_{i, k} (u) = \frac{u - u _{i}}{u _{i + k - 1} - u _{i}} N_{i, k - 1} (u) + \frac{u _{i + k} - u}{u _{i + k} - u _{i + 1}} N_{i + 1, k - 1} (u)$

3.2 主要性质

局部性（Local support）：每个控制点只影响相邻 $k$ 个区间
高阶连续性：可实现 $C^{k - 1}$ 连续
可调节点向量（Knot vector）：灵活控制曲线形状

3.3 de Boor 算法

算法流程

给定节点向量 $U$ ，控制点 $P_{0}, ..., P_{n}$ ，阶数 $k$ ，参数 $u$
找到 $u$ 所在区间 $[u_{j}, u_{j + 1})$
递归线性插值，最终得到 $r (u)$

例题

已知三次 B-spline 曲线，控制点 $P_{0} = (0, 0)$ , $P_{1} = (1, 2)$ , $P_{2} = (3, 3)$ , $P_{3} = (4, 0)$ ，节点向量 $[0, 0, 0, 0, 1, 2, 2, 2, 2]$ ，求 $u = 1$ 处的曲线点。

解答：

节点向量为 $[0, 0, 0, 0, 1, 2, 2, 2, 2]$ ，阶数 $k = 4$ （因为节点向量长度为 $n + k + 1 = 4 + 4 + 1 = 9$ ）。
参数 $u = 1$ 落在区间 $[u_{3}, u_{4}) = [0, 1)$ ，但由于节点向量重复，实际有效区间需要考虑 $u = 1$ 属于 $[u_{4}, u_{5}) = [1, 2)$ ，因此 $j = 4$ 。
影响 $u = 1$ 的控制点为 $P_{1}, P_{2}, P_{3}$ （因为 $j - k + 1 = 4 - 4 + 1 = 1$ 到 $j = 4$ ，即控制点索引 $1$ 到 $3$ ）。
计算基函数 $N_{i, 4} (1)$ ：
- 先计算 $N_{i, 1} (1)$ ：
  - $N_{1, 1} (1) = 0$ （因为 $u = 1$ 不在 $[u_{1}, u_{2}) = [0, 0)$ ）
  - $N_{2, 1} (1) = 0$ （因为 $u = 1$ 不在 $[u_{2}, u_{3}) = [0, 0)$ ）
  - $N_{3, 1} (1) = 0$ （因为 $u = 1$ 不在 $[u_{3}, u_{4}) = [0, 1)$ ）
  - $N_{4, 1} (1) = 1$ （因为 $u = 1$ 在 $[u_{4}, u_{5}) = [1, 2)$ ）
  - $N_{5, 1} (1) = 0$ （因为 $u = 1$ 不在 $[u_{5}, u_{6}) = [2, 2)$ ）
- 递归计算 $N_{i, 2} (1)$ ：
  - $N_{2, 2} (1) = \frac{1 - u _{2}}{u _{3} - u _{2}} N_{2, 1} (1) + \frac{u _{4} - 1}{u _{4} - u _{3}} N_{3, 1} (1) = \frac{1 - 0}{0 - 0} * 0 + \frac{1 - 1}{1 - 0} * 0 = 0$
  - $N_{3, 2} (1) = \frac{1 - u _{3}}{u _{4} - u _{3}} N_{3, 1} (1) + \frac{u _{5} - 1}{u _{5} - u _{4}} N_{4, 1} (1) = \frac{1 - 0}{1 - 0} * 0 + \frac{2 - 1}{2 - 1} * 1 = 1$
  - $N_{4, 2} (1) = \frac{1 - u _{4}}{u _{5} - u _{4}} N_{4, 1} (1) + \frac{u _{6} - 1}{u _{6} - u _{5}} N_{5, 1} (1) = \frac{1 - 1}{2 - 1} * 1 + \frac{2 - 1}{2 - 2} * 0 = 0$
- 递归计算 $N_{i, 3} (1)$ ：
  - $N_{2, 3} (1) = \frac{1 - u _{2}}{u _{4} - u _{2}} N_{2, 2} (1) + \frac{u _{5} - 1}{u _{5} - u _{3}} N_{3, 2} (1) = \frac{1 - 0}{1 - 0} * 0 + \frac{2 - 1}{2 - 0} * 1 = 0.5$
  - $N_{3, 3} (1) = \frac{1 - u _{3}}{u _{5} - u _{3}} N_{3, 2} (1) + \frac{u _{6} - 1}{u _{6} - u _{4}} N_{4, 2} (1) = \frac{1 - 0}{2 - 0} * 1 + \frac{2 - 1}{2 - 1} * 0 = 0.5$
  - $N_{4, 3} (1) = 0$ （类似计算）
- 递归计算 $N_{i, 4} (1)$ ：
  - $N_{1, 4} (1) = \frac{1 - u _{1}}{u _{4} - u _{1}} N_{1, 3} (1) + \frac{u _{5} - 1}{u _{5} - u _{2}} N_{2, 3} (1) = \frac{1 - 0}{1 - 0} * 0 + \frac{2 - 1}{2 - 0} * 0.5 = 0.25$
  - $N_{2, 4} (1) = \frac{1 - u _{2}}{u _{5} - u _{2}} N_{2, 3} (1) + \frac{u _{6} - 1}{u _{6} - u _{3}} N_{3, 3} (1) = \frac{1 - 0}{2 - 0} * 0.5 + \frac{2 - 1}{2 - 0} * 0.5 = 0.5$
  - $N_{3, 4} (1) = \frac{1 - u _{3}}{u _{6} - u _{3}} N_{3, 3} (1) + \frac{u _{7} - 1}{u _{7} - u _{4}} N_{4, 3} (1) = \frac{1 - 0}{2 - 0} * 0.5 + \frac{2 - 1}{2 - 1} * 0 = 0.25$
最终曲线点 $r (1) = N_{1, 4} (1) P_{1} + N_{2, 4} (1) P_{2} + N_{3, 4} (1) P_{3} = 0.25 * (1, 2) + 0.5 * (3, 3) + 0.25 * (4, 0) = (0.25 * 1 + 0.5 * 3 + 0.25 * 4, 0.25 * 2 + 0.5 * 3 + 0.25 * 0) = (2.75, 2.0)$
答： $r (1) = (2.75, 2.0)$

Python代码示例 (手动计算基函数)：

import numpy as np

# 例题数据
control_points_bspline = np.array([[0, 0], [1, 2], [3, 3], [4, 0]])
knots_bspline = np.array([0,0,0,0,1,2,2,2,2])
u_bspline = 1

# 按照文档手动计算的基函数值
N14_at_1 = 0.25
N24_at_1 = 0.5
N34_at_1 = 0.25

# 影响 u=1 的控制点是 P1, P2, P3 (索引1, 2, 3)
point_on_bspline = N14_at_1 * control_points_bspline[1] + \
                   N24_at_1 * control_points_bspline[2] + \
                   N34_at_1 * control_points_bspline[3]
                   
print(f"例题 (B-spline 手动基函数): r({u_bspline}) = {point_on_bspline}")

# 验证 (可选)
# expected_point = np.array([2.75, 2.0])
# assert np.allclose(point_on_bspline, expected_point)

# 注意: 实际应用中会使用递归函数计算基函数或直接使用 scipy.interpolate
# from scipy.interpolate import BSpline
# k = 3 # B样条阶数 (degree), k=order-1. Order k=4 for cubic.
# tck = (knots_bspline, control_points_bspline.T, k)
# point_scipy = BSpline(*tck)(u_bspline)
# print(f"例题 (B-spline Scipy): r({u_bspline}) = {point_scipy}")

新增例题：计算给定控制点 $P_{0} = (0, 0)$ ， $P_{1} = (1, 1)$ ， $P_{2} = (2, 0)$ 的二次Bezier曲线在 $t = 0.5$ 处的点

解答步骤：

使用de Casteljau算法：
- 第一层： $P_{0}^{(1)} = 0.5 \times P_{0} + 0.5 \times P_{1} = (0.5, 0.5)$
- $P_{1}^{(1)} = 0.5 \times P_{1} + 0.5 \times P_{2} = (1.5, 0.5)$
第二层： $P_{0}^{(2)} = 0.5 \times P_{0}^{(1)} + 0.5 \times P_{1}^{(1)} = (1, 0.5)$
结果：曲线点为 $(1, 0.5)$ 。

Python代码示例 (使用上面的 de_casteljau 函数)：

# 新增例题数据
P_bezier_add = [[0, 0], [1, 1], [2, 0]]
t_bezier_add = 0.5

# 计算
point_on_curve_add = de_casteljau(P_bezier_add, t_bezier_add)
print(f"新增例题 (Bezier de Casteljau): P({t_bezier_add}) = {point_on_curve_add}")

# 验证 (可选)
# expected_point = np.array([1, 0.5])
# assert np.allclose(point_on_curve_add, expected_point)

4. 插值方法（Interpolation）

4.1 Lagrange 插值

公式：
$L (x) = i = 0 \sum n y_{i} j \neq = i \prod \frac{x - x _{j}}{x _{i} - x _{j}}$
优缺点：公式显式，适合手算，但高阶时数值不稳定。

例题

已知点 $(0, 1)$ , $(1, 2)$ , $(2, 0)$ ，用 Lagrange 插值求 $x = 1.5$ 处的函数值。

解答：

基函数：
- $l_{0} (x) = \frac{( x - 1 ) ( x - 2 )}{( 0 - 1 ) ( 0 - 2 )} = \frac{( x - 1 ) ( x - 2 )}{2}$
- $l_{1} (x) = \frac{( x - 0 ) ( x - 2 )}{( 1 - 0 ) ( 1 - 2 )} = \frac{x ( x - 2 )}{- 1} = - x (x - 2)$
- $l_{2} (x) = \frac{( x - 0 ) ( x - 1 )}{( 2 - 0 ) ( 2 - 1 )} = \frac{x ( x - 1 )}{2}$
插值多项式： $L (x) = 1 * l_{0} (x) + 2 * l_{1} (x) + 0 * l_{2} (x) = \frac{( x - 1 ) ( x - 2 )}{2} - 2 x (x - 2)$
代入 $x = 1.5$ ：
- $l_{0} (1.5) = \frac{( 1.5 - 1 ) ( 1.5 - 2 )}{2} = \frac{0.5 * ( - 0.5 )}{2} = - 0.125$
- $l_{1} (1.5) = - 1.5 * (1.5 - 2) = - 1.5 * (- 0.5) = 0.75$
- $L (1.5) = 1 * (- 0.125) + 2 * (0.75) = - 0.125 + 1.5 = 1.375$
答： $L (1.5) = 1.375$

Python代码示例 (Lagrange 插值)：

import numpy as np
from scipy import interpolate # 用于验证

def lagrange_basis(x_data, i, x):
    """计算第 i 个 Lagrange 基函数 l_i(x)"""
    n = len(x_data)
    li = 1.0
    for j in range(n):
        if i != j:
            li *= (x - x_data[j]) / (x_data[i] - x_data[j])
    return li

def lagrange_interpolation(x_data, y_data, x_interp):
    """使用 Lagrange 插值计算 x_interp 处的函数值"""
    n = len(x_data)
    result = 0.0
    for i in range(n):
        result += y_data[i] * lagrange_basis(x_data, i, x_interp)
    return result

# 例题数据
x_lagrange = np.array([0, 1, 2])
y_lagrange = np.array([1, 2, 0])
x_interp_lagrange = 1.5

# 计算
interp_value_lagrange = lagrange_interpolation(x_lagrange, y_lagrange, x_interp_lagrange)
print(f"例题 (Lagrange 插值): L({x_interp_lagrange}) = {interp_value_lagrange}")

# 验证 (可选)
# expected_value = 1.375
# assert np.isclose(interp_value_lagrange, expected_value)
# lagrange_poly = interpolate.lagrange(x_lagrange, y_lagrange)
# scipy_value = lagrange_poly(x_interp_lagrange)
# assert np.isclose(scipy_value, expected_value)

4.2 Newton 插值

公式：
$N (x) = a_{0} + a_{1} (x - x_{0}) + a_{2} (x - x_{0}) (x - x_{1}) + \dots$ 其中 $a_{i}$ 为差商（divided differences）

例题

已知点 $(0, 1)$ , $(1, 2)$ , $(2, 0)$ ，用 Newton 插值求 $x = 1.5$ 处的函数值。

解答：

差商表：
- $f [0] = 1$
- $f [1] = 2$
- $f [2] = 0$
- $f [0, 1] = \frac{f [ 1 ] - f [ 0 ]}{1 - 0} = \frac{2 - 1}{1} = 1$
- $f [1, 2] = \frac{f [ 2 ] - f [ 1 ]}{2 - 1} = \frac{0 - 2}{1} = - 2$
- $f [0, 1, 2] = \frac{f [ 1 , 2 ] - f [ 0 , 1 ]}{2 - 0} = \frac{- 2 - 1}{2} = - 1.5$
插值多项式： $N (x) = 1 + 1 * (x - 0) + (- 1.5) * (x - 0) (x - 1)$
代入 $x = 1.5$ ：
- $N (1.5) = 1 + 1 * (1.5) + (- 1.5) * (1.5) * (1.5 - 1) = 1 + 1.5 + (- 1.5) * (1.5) * (0.5) = 1 + 1.5 - 1.125 = 1.375$
答： $N (1.5) = 1.375$

Python代码示例 (Newton 插值)：

import numpy as np

def divided_differences(x_data, y_data):
    """计算差商表，返回 Newton 多项式系数"""
    n = len(x_data)
    coef = np.zeros([n, n])
    coef[:, 0] = y_data # 第一列是 y 值
    
    for j in range(1, n): # 列
        for i in range(n - j): # 行
            coef[i, j] = (coef[i + 1, j - 1] - coef[i, j - 1]) / (x_data[i + j] - x_data[i])
            
    return coef[0, :] # 返回第一行作为系数 a0, a1, a2, ...

def newton_polynomial(coefficients, x_data, x):
    """根据差商系数计算 Newton 插值多项式在 x 处的值"""
    n = len(coefficients)
    result = coefficients[0]
    term = 1.0
    
    for i in range(1, n):
        term *= (x - x_data[i - 1]) # 计算 (x-x0), (x-x0)(x-x1), ...
        result += coefficients[i] * term
        
    return result

# 例题数据
x_newton = np.array([0, 1, 2])
y_newton = np.array([1, 2, 0])
x_interp_newton = 1.5

# 计算
coeffs_newton = divided_differences(x_newton, y_newton)
interp_value_newton = newton_polynomial(coeffs_newton, x_newton, x_interp_newton)
print(f"例题 (Newton 插值) 差商系数: {coeffs_newton}")
print(f"例题 (Newton 插值): N({x_interp_newton}) = {interp_value_newton}")

# 验证 (可选)
# expected_value = 1.375
# assert np.isclose(interp_value_newton, expected_value)

4.3 样条插值（Spline Interpolation）

三次样条：分段三次多项式，保证 $C^{2}$ 连续，常用于平滑曲线拟合。

例题

已知点 $(0, 1)$ , $(1, 2)$ , $(2, 0)$ ，构造自然三次样条插值（边界条件为二阶导数为0）。

解答：

设区间 $[0, 1]$ 和 $[1, 2]$ 上的三次多项式为：
- $S_{1} (x) = a_{1} + b_{1} (x - 0) + c_{1} (x - 0)^{2} + d_{1} (x - 0)^{3}$
- $S_{2} (x) = a_{2} + b_{2} (x - 1) + c_{2} (x - 1)^{2} + d_{2} (x - 1)^{3}$
条件：
1. 插值条件： $S_{1} (0) = 1$ , $S_{1} (1) = 2$ , $S_{2} (1) = 2$ , $S_{2} (2) = 0$
2. 一阶导数连续： $S_{1}^{'} (1) = S_{2}^{'} (1)$
3. 二阶导数连续： $S_{1}^{''} (1) = S_{2}^{''} (1)$
4. 自然边界条件： $S_{1}^{''} (0) = 0$ , $S_{2}^{''} (2) = 0$
解方程组（详细计算略）：
- $S_{1} (x) = 1 + 1.5 x - 0.5 x^{3}$
- $S_{2} (x) = 2 - 0.5 (x - 1) - 1.5 (x - 1)^{2} + 0.5 (x - 1)^{3}$
答：自然三次样条插值多项式如上。

Python代码示例 (使用 Scipy)：

import numpy as np
from scipy.interpolate import CubicSpline

# 例题数据
x_spline = np.array([0, 1, 2])
y_spline = np.array([1, 2, 0])

# 计算自然三次样条 (bc_type='natural')
cs = CubicSpline(x_spline, y_spline, bc_type='natural')

# 打印样条系数 (可选)
# print("样条系数 (分段):")
# for i in range(len(cs.c[0, :])):
#     print(f" 区间 {i}: {cs.c[:, i]}")

# 在 x=1.5 处插值 (虽然例题只要求构造)
x_interp_spline = 1.5
interp_value_spline = cs(x_interp_spline)
print(f"例题 (自然三次样条插值): S({x_interp_spline}) = {interp_value_spline}")

# 验证 x=0, 1, 2 处的插值结果
print(f"验证: S(0)={cs(0)}, S(1)={cs(1)}, S(2)={cs(2)}")
assert np.isclose(cs(0), 1)
assert np.isclose(cs(1), 2)
assert np.isclose(cs(2), 0)

# 验证边界条件 S''(0)=0, S''(2)=0
print(f"验证: S''(0)={cs(0, 2)}, S''(2)={cs(2, 2)}")
assert np.isclose(cs(0, 2), 0)
assert np.isclose(cs(2, 2), 0)

5. 常见考点与易错点

Bezier 曲线的端点、切线、凸包等性质
de Casteljau 算法步骤与递归思想
B-spline 节点向量与控制点影响范围
插值多项式的数值稳定性与误差
英文术语拼写与公式记忆

6. 英文关键词与表达

Bezier curve, control point, Bernstein polynomial, de Casteljau algorithm, subdivision, convex hull, affine invariance
B-spline, knot vector, de Boor point, local support, continuity, degree, order
interpolation, Lagrange interpolation, Newton divided difference, spline, cubic spline

如需更详细例题推导或某一算法的代码实现，可随时补充！

CE7453 Numerical Algorithms 期末高分超详细攻略（第四章：Numerical Differentiation & Integration）

本章内容：数值微分、数值积分
适用对象：零基础/考前冲刺/快速查漏补缺
内容特色：详细原理、公式推导、算法流程、例题全解、英文关键词

1. 基本概念与背景

数值微分（Numerical Differentiation）：用离散点近似计算函数的导数，适用于函数表达式复杂或仅有数据点的情况。
数值积分（Numerical Integration）：用有限和近似计算定积分，适用于无法解析积分或仅有数据点的情况。
实际应用：工程仿真、信号处理、物理建模、数据分析等。

2. 数值微分（Numerical Differentiation）

2.1 差分公式（Finite Difference Formulas）

2.1.1 2点前向差分（2-point forward difference）

公式：
$f^{'} (x) \approx \frac{f ( x + h ) - f ( x )}{h}$
误差： $O (h)$ ，截断误差较大

2.1.2 3点中心差分（3-point centered difference）

公式：
$f^{'} (x) \approx \frac{f ( x + h ) - f ( x - h )}{2 h}$
误差： $O (h^{2})$ ，更精确

2.1.3 5点中心差分（5-point centered difference）

公式：
$f^{'} (x) \approx \frac{f ( x - 2 h ) - 8 f ( x - h ) + 8 f ( x + h ) - f ( x + 2 h )}{12 h}$
误差： $O (h^{4})$ ，高精度

2.1.4 Richardson 外推（Richardson Extrapolation）

思想：用不同步长的近似结果组合，消除主误差项，提高精度。

2.2 算法流程（以3点中心差分为例）

选定步长 $h$ （过大误差大，过小舍入误差大）
计算 $f (x + h)$ 和 $f (x - h)$
按公式求导数近似值
可用不同 $h$ 验证结果稳定性

2.3 典型例题

例题1：已知 $f (x) = e^{2 x}$ ，用 $h = 0.1$ ，求 $x = 1$ 处的导数近似值（用前向差分法）

解答：

$f (1) = e^{2 \cdot 1} = e^{2} \approx 7.3891$
$f (1 + 0.1) = e^{2 \cdot 1.1} = e^{2.2} \approx 9.0250$
$f^{'} (1) \approx \frac{9.0250 - 7.3891}{0.1} = \frac{1.6359}{0.1} \approx 16.3596$
精确值 $f^{'} (1) = 2 e^{2} \approx 14.7782$ ，误差约为 $1.5814$

例题2：已知 $f (x) = e^{2 x}$ ，用 $h = 0.1$ ，求 $x = 1$ 处的导数近似值（用中心差分法）

解答：

$f (1 + 0.1) = e^{2 \cdot 1.1} = e^{2.2} \approx 9.0250$
$f (1 - 0.1) = e^{2 \cdot 0.9} = e^{1.8} \approx 6.0496$
$f^{'} (1) \approx \frac{9.0250 - 6.0496}{2 \cdot 0.1} = \frac{2.9754}{0.2} \approx 14.8768$
精确值 $f^{'} (1) = 2 e^{2} \approx 14.7782$ ，误差约为 $0.0986$

例题3：已知 $f (x) = x^{3} + 2 x^{2} + 1$ ，在 $x = 2$ 处用 $h = 0.2$ ，分别用2点前向差分、3点中心差分和5点中心差分计算导数近似值，并与精确值比较。

解答：

精确导数： $f^{'} (x) = 3 x^{2} + 4 x$ ，在 $x = 2$ 处， $f^{'} (2) = 3 * (2)^{2} + 4 * 2 = 12 + 8 = 20$
2点前向差分：
- $f (2) = 2^{3} + 2 * 2^{2} + 1 = 8 + 8 + 1 = 17$
- $f (2 + 0.2) = f (2.2) = (2.2)^{3} + 2 * (2.2)^{2} + 1 = 10.648 + 9.68 + 1 = 21.328$
- $f^{'} (2) \approx \frac{21.328 - 17}{0.2} = \frac{4.328}{0.2} = 21.64$
- 误差： $∣21.64 - 20∣ = 1.64$
3点中心差分：
- $f (2 - 0.2) = f (1.8) = (1.8)^{3} + 2 * (1.8)^{2} + 1 = 5.832 + 6.48 + 1 = 13.312$
- $f (2 + 0.2) = 21.328$ （如上）
- $f^{'} (2) \approx \frac{21.328 - 13.312}{2 * 0.2} = \frac{8.016}{0.4} = 20.04$
- 误差： $∣20.04 - 20∣ = 0.04$
5点中心差分：
- $f (2 - 0.4) = f (1.6) = (1.6)^{3} + 2 * (1.6)^{2} + 1 = 4.096 + 5.12 + 1 = 10.216$
- $f (2 - 0.2) = 13.312$ （如上）
- $f (2 + 0.2) = 21.328$ （如上）
- $f (2 + 0.4) = f (2.4) = (2.4)^{3} + 2 * (2.4)^{2} + 1 = 13.824 + 11.52 + 1 = 26.344$
- $f^{'} (2) \approx \frac{10.216 - 8 * 13.312 + 8 * 21.328 - 26.344}{12 * 0.2}$
- $= \frac{10.216 - 106.496 + 170.624 - 26.344}{2.4} = \frac{47.996}{2.4} = 19.9983$
- 误差： $∣19.9983 - 20∣ = 0.0017$
比较：5点中心差分精度最高，误差最小；2点前向差分误差最大。

2.4 实际应用与注意事项

应用场景：数值微分常用于速度、加速度计算（如运动学分析），以及优化算法中的梯度计算。
注意事项：
- 步长 $h$ 选择需平衡截断误差和舍入误差，建议多次尝试不同 $h$ 值。
- 对于噪声数据，数值微分会放大噪声，需先平滑处理数据。

3. 数值积分（Numerical Integration）

3.1 Newton-Cotes 公式

3.1.1 梯形法（Trapezoidal Rule）

公式：
$\int_{a}^{b} f (x) d x \approx \frac{h}{2} [f (a) + f (b)]$ 其中 $h = b - a$
复合梯形法：将区间分为 $n$ 段，累加每段梯形面积 $\int_{a}^{b} f (x) d x \approx \frac{h}{2} [f (a) + 2 k = 1 \sum n - 1 f (a + kh) + f (b)]$

3.1.2 Simpson 法（Simpson's Rule）

公式：
$\int_{a}^{b} f (x) d x \approx \frac{h}{3} [f (a) + 4 f ((a + b) /2) + f (b)]$ $h = (b - a) /2$
复合Simpson法：区间分偶数段，累加每段Simpson近似 $\int_{a}^{b} f (x) d x \approx \frac{h}{3} [f (a) + 4 k = 1, 3, ... \sum n - 1 f (a + kh) + 2 k = 2, 4, ... \sum n - 2 f (a + kh) + f (b)]$

3.1.3 Romberg 积分（Romberg Integration）

思想：用梯形法结果递推外推，提高精度。基于 Richardson 外推，通过多次梯形法计算，逐步消除误差项。

3.1.4 高斯求积（Gaussian Quadrature）

思想：选取最优节点和权重，使多项式积分精度最高。对于 $n$ 点高斯求积，对最高 $2 n - 1$ 次多项式精确。

3.2 算法流程（以复合梯形法为例）

将 $[a, b]$ 分为 $n$ 段，步长 $h = (b - a) / n$
计算端点 $f (a), f (b)$
计算中间点 $f (a + kh)$ ， $k = 1, 2, ..., n - 1$
按公式累加求和

3.3 典型例题

例题1：用复合梯形法计算 $\int_{0}^{1} e^{x} d x$ ，分4段

解答：

$h = 0.25$
$x_{0} = 0, x_{1} = 0.25, x_{2} = 0.5, x_{3} = 0.75, x_{4} = 1$
$f (x_{0}) = 1, f (x_{1}) = 1.2840, f (x_{2}) = 1.6487, f (x_{3}) = 2.1170, f (x_{4}) = 2.7183$
$I \approx \frac{0.25}{2} [1 + 2 (1.2840 + 1.6487 + 2.1170) + 2.7183] = 1.7183$
精确值 $e - 1 = 1.7183$ ，误差约为 $0.0000$

例题2：用Simpson法计算积分 $\int_{0}^{1} x^{2} d x$ ，分2段

解答步骤：

$h = 0.5$ ，点： $x_{0} = 0, x_{1} = 0.5, x_{2} = 1$
$f (x_{0}) = 0, f (x_{1}) = 0.25, f (x_{2}) = 1$
$I \approx \frac{0.5}{3} [0 + 4 * 0.25 + 1] = \frac{0.5}{3} [0 + 1 + 1] = \frac{0.5}{3} * 2 = \frac{1}{3}$
精确值 $\frac{1}{3}$ ，误差为0。

例题3：用Romberg积分计算 $\int_{0}^{1} \frac{4}{1 + x ^{2}} d x$ ，取初始 $n = 1$ ，进行两次外推。

解答步骤：

该积分精确值为 $π \approx 3.1415926535$ ，用于验证误差。
第一步：梯形法，不同分割：
- $n = 1$ ， $h = 1$ ， $T_{1} (1) = \frac{1}{2} [f (0) + f (1)] = \frac{1}{2} [4 + 2] = 3$
- $n = 2$ ， $h = 0.5$ ， $T_{1} (2) = \frac{0.5}{2} [f (0) + 2 f (0.5) + f (1)] = \frac{0.5}{2} [4 + 2 * 3.2 + 2] = 3.2$
- $n = 4$ ， $h = 0.25$ ， $T_{1} (4) = \frac{0.25}{2} [f (0) + 2 (f (0.25) + f (0.5) + f (0.75)) + f (1)]$
  - $f (0.25) = \frac{4}{1 + 0.0625} = 3.7647$ , $f (0.5) = 3.2$ , $f (0.75) = \frac{4}{1 + 0.5625} = 2.56$
  - $T_{1} (4) = \frac{0.25}{2} [4 + 2 * (3.7647 + 3.2 + 2.56) + 2] = 3.1466$
第二步：Romberg外推：
- 第一次外推： $T_{2} (2) = \frac{4 * T _{1} ( 2 ) - T _{1} ( 1 )}{4 - 1} = \frac{4 * 3.2 - 3}{3} = 3.2667$
- 第二次外推： $T_{2} (4) = \frac{4 * T _{1} ( 4 ) - T _{1} ( 2 )}{4 - 1} = \frac{4 * 3.1466 - 3.2}{3} = 3.1287$
- 继续外推： $T_{3} (4) = \frac{4 * T _{2} ( 4 ) - T _{2} ( 2 )}{4 - 1} = \frac{4 * 3.1287 - 3.2667}{3} = 3.1493$
结果： $T_{3} (4) = 3.1493$ ，与 $π$ 较接近，误差约为 $0.0077$ 。

例题4：用三点高斯求积法计算 $\int_{- 1}^{1} 1 - x^{2} d x$

解答步骤：

三点高斯求积法使用的节点： $x_{1} = - 3/5$ ， $x_{2} = 0$ ， $x_{3} = 3/5$
对应的权重： $w_{1} = 5/9$ ， $w_{2} = 8/9$ ， $w_{3} = 5/9$
函数值计算：
- $f (x_{1}) = f (- 3/5) = 1 - (3/5)^{2} = 1 - 3/5 = 2/5 \approx 0.6325$
- $f (x_{2}) = f (0) = 1 - 0^{2} = 1$
- $f (x_{3}) = f (3/5) = 1 - (3/5)^{2} = 1 - 3/5 = 2/5 \approx 0.6325$
求积结果： $I \approx w_{1} f (x_{1}) + w_{2} f (x_{2}) + w_{3} f (x_{3})$
- $I \approx (5/9) \cdot 0.6325 + (8/9) \cdot 1 + (5/9) \cdot 0.6325$
- $I \approx (5/9) \cdot 0.6325 \cdot 2 + (8/9) \cdot 1$
- $I \approx 0.7028 + 0.8889 \approx 1.5916$
精确值为 $π /2 \approx 1.5708$ ，误差约为 $0.0208$ （约1.3%）

3.4 实际应用与注意事项

应用场景：数值积分用于计算面积、体积、概率密度积分、物理量累积（如功、能量）。
注意事项：
- 增加分割数 $n$ 可提高精度，但计算量增加。
- 对于周期性函数或振荡函数，高斯求积可能更有效。
- Romberg 积分适合快速收敛，但对函数平滑性要求较高。

4. 常见考点与易错点

步长 $h$ 选取不当导致误差大
舍入误差与截断误差权衡
复合公式累加时端点权重
高斯求积节点和权重记忆
英文术语拼写与公式记忆

5. 英文关键词与表达

numerical differentiation, finite difference, truncation error, round-off error, Richardson extrapolation, step size
numerical integration, Newton-Cotes, trapezoidal rule, Simpson's rule, Romberg integration, Gaussian quadrature, composite rule, node, weight

6. 代码实现示例

下面提供各种数值微分和数值积分方法的Python实现示例。

6.1 数值微分代码实现

import numpy as np
import matplotlib.pyplot as plt

def forward_difference(f, x, h=0.01):
    """
    使用前向差分法计算导数
    
    参数:
        f: 函数
        x: 计算导数的点
        h: 步长
    
    返回:
        导数近似值
    """
    return (f(x + h) - f(x)) / h

def central_difference(f, x, h=0.01):
    """
    使用中心差分法计算导数
    
    参数:
        f: 函数
        x: 计算导数的点
        h: 步长
    
    返回:
        导数近似值
    """
    return (f(x + h) - f(x - h)) / (2 * h)

def five_point_difference(f, x, h=0.01):
    """
    使用五点中心差分法计算导数
    
    参数:
        f: 函数
        x: 计算导数的点
        h: 步长
    
    返回:
        导数近似值
    """
    return (f(x - 2*h) - 8*f(x - h) + 8*f(x + h) - f(x + 2*h)) / (12 * h)

def richardson_extrapolation(f, x, h=0.1, k=2):
    """
    使用Richardson外推法提高导数计算精度
    
    参数:
        f: 函数
        x: 计算导数的点
        h: 初始步长
        k: 外推次数
    
    返回:
        改进后的导数近似值
    """
    # 计算中心差分
    D1 = central_difference(f, x, h)
    D2 = central_difference(f, x, h/2)
    
    # 外推公式: D = D2 + (D2 - D1)/(4^k - 1)
    return D2 + (D2 - D1) / (4**k - 1)

# 测试示例
if __name__ == "__main__":
    # 定义测试函数和其精确导数
    f = lambda x: np.exp(2*x)
    df_exact = lambda x: 2 * np.exp(2*x)
    
    # 测试点
    x0 = 1.0
    
    # 计算导数并与精确值比较
    h = 0.1
    fd = forward_difference(f, x0, h)
    cd = central_difference(f, x0, h)
    fpd = five_point_difference(f, x0, h)
    rd = richardson_extrapolation(f, x0, h)
    exact = df_exact(x0)
    
    print(f"函数 f(x) = e^(2x) 在 x = {x0} 处的导数:")
    print(f"精确值: {exact:.6f}")
    print(f"前向差分 (h={h}): {fd:.6f}, 相对误差: {abs(fd-exact)/exact:.6f}")
    print(f"中心差分 (h={h}): {cd:.6f}, 相对误差: {abs(cd-exact)/exact:.6f}")
    print(f"五点差分 (h={h}): {fpd:.6f}, 相对误差: {abs(fpd-exact)/exact:.6f}")
    print(f"Richardson外推: {rd:.6f}, 相对误差: {abs(rd-exact)/exact:.6f}")
    
    # 绘制不同步长下的误差
    h_values = np.logspace(-8, -1, 20)
    fd_errors = [abs(forward_difference(f, x0, h) - exact)/exact for h in h_values]
    cd_errors = [abs(central_difference(f, x0, h) - exact)/exact for h in h_values]
    fpd_errors = [abs(five_point_difference(f, x0, h) - exact)/exact for h in h_values]
    
    plt.figure(figsize=(10, 6))
    plt.loglog(h_values, fd_errors, 'o-', label='前向差分')
    plt.loglog(h_values, cd_errors, 's-', label='中心差分')
    plt.loglog(h_values, fpd_errors, '^-', label='五点差分')
    plt.grid(True)
    plt.xlabel('步长 h')
    plt.ylabel('相对误差')
    plt.title('不同数值微分方法在不同步长下的相对误差')
    plt.legend()
    plt.show()

6.2 数值积分代码实现

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

def trapezoidal_rule(f, a, b, n=100):
    """
    复合梯形法则计算定积分
    
    参数:
        f: 被积函数
        a, b: 积分区间
        n: 分段数
    
    返回:
        积分近似值
    """
    h = (b - a) / n
    x = np.linspace(a, b, n+1)
    y = f(x)
    
    # 复合梯形法则公式
    integral = (h/2) * (y[0] + 2*np.sum(y[1:-1]) + y[-1])
    return integral

def simpson_rule(f, a, b, n=100):
    """
    复合Simpson法则计算定积分
    
    参数:
        f: 被积函数
        a, b: 积分区间
        n: 分段数 (必须为偶数)
    
    返回:
        积分近似值
    """
    if n % 2 != 0:
        n += 1  # 确保n为偶数
    
    h = (b - a) / n
    x = np.linspace(a, b, n+1)
    y = f(x)
    
    # 复合Simpson法则公式
    integral = (h/3) * (y[0] + 4*np.sum(y[1:-1:2]) + 2*np.sum(y[2:-1:2]) + y[-1])
    return integral

def romberg_integration(f, a, b, max_iter=5):
    """
    Romberg积分法
    
    参数:
        f: 被积函数
        a, b: 积分区间
        max_iter: 最大迭代次数
        
    返回:
        积分近似值
    """
    R = np.zeros((max_iter, max_iter))
    
    # 初始计算 - 使用复合梯形法则
    h = b - a
    R[0, 0] = (h/2) * (f(a) + f(b))
    
    for i in range(1, max_iter):
        # 计算下一级梯形法值
        h = h / 2
        n = 2**i
        sum_f = 0
        for k in range(1, n, 2):
            sum_f += f(a + k*h)
        R[i, 0] = R[i-1, 0]/2 + h*sum_f
        
        # Richardson外推
        for j in range(1, i+1):
            R[i, j] = R[i, j-1] + (R[i, j-1] - R[i-1, j-1]) / (4**j - 1)
    
    return R[max_iter-1, max_iter-1]

def gaussian_quadrature(f, a, b, n=5):
    """
    高斯求积法计算定积分
    
    参数:
        f: 被积函数
        a, b: 积分区间
        n: 高斯点数 (最多5点)
        
    返回:
        积分近似值
    """
    # 高斯点和权重 (标准区间 [-1, 1])
    if n == 1:
        x_points = np.array([0])
        weights = np.array([2])
    elif n == 2:
        x_points = np.array([-1/np.sqrt(3), 1/np.sqrt(3)])
        weights = np.array([1, 1])
    elif n == 3:
        x_points = np.array([-np.sqrt(3/5), 0, np.sqrt(3/5)])
        weights = np.array([5/9, 8/9, 5/9])
    elif n == 4:
        x0 = np.sqrt((3-2*np.sqrt(6/5))/7)
        x1 = np.sqrt((3+2*np.sqrt(6/5))/7)
        x_points = np.array([-x1, -x0, x0, x1])
        w0 = (18+np.sqrt(30))/36
        w1 = (18-np.sqrt(30))/36
        weights = np.array([w1, w0, w0, w1])
    elif n == 5:
        x0 = 0
        x1 = np.sqrt(5-2*np.sqrt(10/7))/3
        x2 = np.sqrt(5+2*np.sqrt(10/7))/3
        x_points = np.array([-x2, -x1, x0, x1, x2])
        w0 = 128/225
        w1 = (322+13*np.sqrt(70))/900
        w2 = (322-13*np.sqrt(70))/900
        weights = np.array([w2, w1, w0, w1, w2])
    else:
        raise ValueError("目前仅支持1到5个高斯点")
    
    # 转换到积分区间 [a, b]
    scaled_f = lambda x: f((b-a)/2 * x + (a+b)/2) * (b-a)/2
    
    # 计算积分
    integral = np.sum(weights * scaled_f(x_points))
    return integral

# 测试示例
if __name__ == "__main__":
    # 测试函数
    f1 = lambda x: x**2
    f2 = lambda x: np.exp(x)
    f3 = lambda x: 4/(1+x**2)  # 积分结果为π
    f4 = lambda x: np.sqrt(1-x**2)  # 半圆面积，积分结果为π/2
    
    # 积分区间
    a1, b1 = 0, 1  # 用于f1, f2
    a3, b3 = 0, 1  # 用于f3
    a4, b4 = -1, 1  # 用于f4
    
    # 精确结果
    exact1 = 1/3  # ∫(0→1) x^2 dx = 1/3
    exact2 = np.exp(1) - 1  # ∫(0→1) e^x dx = e - 1
    exact3 = np.pi  # ∫(0→1) 4/(1+x^2) dx = π
    exact4 = np.pi/2  # ∫(-1→1) sqrt(1-x^2) dx = π/2
    
    # 不同方法计算结果
    n_intervals = 4  # 分段数
    
    print("函数 f(x) = x^2 在区间 [0, 1] 上的积分:")
    trap1 = trapezoidal_rule(f1, a1, b1, n_intervals)
    simp1 = simpson_rule(f1, a1, b1, n_intervals)
    romb1 = romberg_integration(f1, a1, b1, 3)
    gauss1 = gaussian_quadrature(f1, a1, b1, 3)
    
    print(f"精确值: {exact1}")
    print(f"梯形法则 (n={n_intervals}): {trap1:.6f}, 相对误差: {abs(trap1-exact1)/exact1:.6f}")
    print(f"Simpson法则 (n={n_intervals}): {simp1:.6f}, 相对误差: {abs(simp1-exact1)/exact1:.6f}")
    print(f"Romberg积分: {romb1:.6f}, 相对误差: {abs(romb1-exact1)/exact1:.6f}")
    print(f"高斯求积法 (3点): {gauss1:.6f}, 相对误差: {abs(gauss1-exact1)/exact1:.6f}")
    
    print("\n函数 f(x) = 4/(1+x^2) 在区间 [0, 1] 上的积分 (结果应接近π):")
    trap3 = trapezoidal_rule(f3, a3, b3, n_intervals)
    simp3 = simpson_rule(f3, a3, b3, n_intervals)
    romb3 = romberg_integration(f3, a3, b3, 3)
    gauss3 = gaussian_quadrature(f3, a3, b3, 3)
    
    print(f"精确值: {exact3}")
    print(f"梯形法则 (n={n_intervals}): {trap3:.6f}, 相对误差: {abs(trap3-exact3)/exact3:.6f}")
    print(f"Simpson法则 (n={n_intervals}): {simp3:.6f}, 相对误差: {abs(simp3-exact3)/exact3:.6f}")
    print(f"Romberg积分: {romb3:.6f}, 相对误差: {abs(romb3-exact3)/exact3:.6f}")
    print(f"高斯求积法 (3点): {gauss3:.6f}, 相对误差: {abs(gauss3-exact3)/exact3:.6f}")
    
    # 绘制不同分段数下的误差比较
    n_values = [2, 4, 8, 16, 32, 64, 128]
    
    trap_errors = [abs(trapezoidal_rule(f2, a1, b1, n) - exact2)/exact2 for n in n_values]
    simp_errors = [abs(simpson_rule(f2, a1, b1, n) - exact2)/exact2 for n in n_values]
    
    plt.figure(figsize=(10, 6))
    plt.loglog(n_values, trap_errors, 'o-', label='梯形法则')
    plt.loglog(n_values, simp_errors, 's-', label='Simpson法则')
    plt.grid(True)
    plt.xlabel('分段数 n')
    plt.ylabel('相对误差')
    plt.title('不同数值积分方法在不同分段数下的相对误差 (f(x) = e^x)')
    plt.legend()
    plt.show()

6.3 应用实例：微分方程数值解法

import numpy as np
import matplotlib.pyplot as plt

def euler_method(f, t0, y0, h, n_steps):
    """
    欧拉法求解常微分方程 dy/dt = f(t, y)
    
    参数:
        f: 函数，表示 dy/dt = f(t, y)
        t0: 初始时间
        y0: 初始值
        h: 步长
        n_steps: 步数
    
    返回:
        t_values: 时间点数组
        y_values: 对应的函数值数组
    """
    t_values = np.zeros(n_steps + 1)
    y_values = np.zeros(n_steps + 1)
    
    t_values[0] = t0
    y_values[0] = y0
    
    for i in range(n_steps):
        t_values[i+1] = t_values[i] + h
        y_values[i+1] = y_values[i] + h * f(t_values[i], y_values[i])
    
    return t_values, y_values

def runge_kutta_4(f, t0, y0, h, n_steps):
    """
    四阶Runge-Kutta方法求解常微分方程 dy/dt = f(t, y)
    
    参数:
        f: 函数，表示 dy/dt = f(t, y)
        t0: 初始时间
        y0: 初始值
        h: 步长
        n_steps: 步数
    
    返回:
        t_values: 时间点数组
        y_values: 对应的函数值数组
    """
    t_values = np.zeros(n_steps + 1)
    y_values = np.zeros(n_steps + 1)
    
    t_values[0] = t0
    y_values[0] = y0
    
    for i in range(n_steps):
        t = t_values[i]
        y = y_values[i]
        
        k1 = f(t, y)
        k2 = f(t + h/2, y + h*k1/2)
        k3 = f(t + h/2, y + h*k2/2)
        k4 = f(t + h, y + h*k3)
        
        t_values[i+1] = t + h
        y_values[i+1] = y + h * (k1 + 2*k2 + 2*k3 + k4) / 6
    
    return t_values, y_values

# 示例：求解微分方程 dy/dt = -y + t (精确解为 y = t - 1 + 2e^(-t))
if __name__ == "__main__":
    # 定义微分方程 dy/dt = f(t, y)
    f = lambda t, y: -y + t
    
    # 定义精确解
    exact_solution = lambda t: t - 1 + 2 * np.exp(-t)
    
    # 初始条件和求解参数
    t0, y0 = 0, 1  # 初始条件
    h = 0.1        # 步长
    n_steps = 50   # 步数
    
    # 使用欧拉法求解
    t_euler, y_euler = euler_method(f, t0, y0, h, n_steps)
    
    # 使用Runge-Kutta法求解
    t_rk4, y_rk4 = runge_kutta_4(f, t0, y0, h, n_steps)
    
    # 计算精确解
    t_exact = np.linspace(t0, t0 + h*n_steps, 200)
    y_exact = exact_solution(t_exact)
    
    # 绘制结果比较
    plt.figure(figsize=(12, 6))
    plt.plot(t_exact, y_exact, 'k-', label='精确解')
    plt.plot(t_euler, y_euler, 'bo-', label='欧拉法')
    plt.plot(t_rk4, y_rk4, 'rs-', label='四阶Runge-Kutta法')
    plt.grid(True)
    plt.xlabel('t')
    plt.ylabel('y')
    plt.title('常微分方程 dy/dt = -y + t 的数值解 (h = ' + str(h) + ')')
    plt.legend()
    
    # 计算并显示误差
    y_exact_at_points = exact_solution(t_euler)
    euler_error = np.abs(y_euler - y_exact_at_points)
    rk4_error = np.abs(y_rk4 - y_exact_at_points)
    
    plt.figure(figsize=(12, 6))
    plt.semilogy(t_euler, euler_error, 'bo-', label='欧拉法误差')
    plt.semilogy(t_euler, rk4_error, 'rs-', label='四阶Runge-Kutta法误差')
    plt.grid(True)
    plt.xlabel('t')
    plt.ylabel('误差 (对数尺度)')
    plt.title('数值解误差比较')
    plt.legend()
    
    plt.show()

如需更详细例题推导或某一算法的代码实现，可随时补充！

CE7453 Numerical Algorithms 期末高分超详细攻略（第五章：Least Squares）

本章内容：最小二乘法、正规方程、QR分解、非线性最小二乘、GPS应用
适用对象：零基础/考前冲刺/快速查漏补缺
内容特色：详细原理、公式推导、算法流程、例题全解、英文关键词

1. 基本概念与背景

最小二乘法（Least Squares）：用于拟合数据、解超定方程组（方程数多于未知数），使误差平方和最小。
实际应用：数据拟合、回归分析、信号处理、GPS定位、机器学习等。

2. 线性最小二乘法

2.1 问题描述

给定 $m$ 个方程 $A x \approx b$ ， $A$ 为 $m \times n$ 矩阵（ $m > n$ ），通常无精确解。
目标：找到 $x$ 使 $∣∣ A x - b ∣ ∣^{2}$ 最小。

2.2 正规方程（Normal Equation）

推导： $x min ∣∣ A x - b ∣ ∣^{2} = (A x - b)^{T} (A x - b)$ 对 $x$ 求导，令导数为0，得 $A^{T} A x = A^{T} b$
解法：解 $n$ 阶方程组 $A^{T} A x = A^{T} b$

2.3 QR分解法（QR Decomposition）

原理：将 $A$ 分解为 $A = QR$ ， $Q$ 为正交矩阵， $R$ 为上三角矩阵。
步骤：
1. $A = QR$
2. $R x = Q^{T} b$ ，回代求解 $x$
优点：数值稳定性高，适合大规模问题

2.4 典型例题

例题1：用最小二乘法拟合直线 $y = a x + b$ ，已知数据点 $(1, 2), (2, 3), (3, 5)$

详细解答：

问题建模：目标是找到参数 $a$ 和 $b$ ，使得直线 $y = a x + b$ 尽可能接近给定的数据点。每个数据点对应一个方程：
- 对于 $(1, 2)$ ： $a \cdot 1 + b = 2$
- 对于 $(2, 3)$ ： $a \cdot 2 + b = 3$
- 对于 $(3, 5)$ ： $a \cdot 3 + b = 5$
矩阵形式：将上述方程组写成矩阵形式 $A x = b$ ： $123111 [a b] = 235$ 这里 $A$ 是 $3 \times 2$ 矩阵， $x = [a, b]^{T}$ 是待求参数向量， $b = [2, 3, 5]^{T}$ 是观测值向量。
正规方程：由于 $A$ 不是方阵，无法直接求逆，我们使用正规方程 $A^{T} A x = A^{T} b$ ：
- 计算 $A^{T}$ ： $A^{T} = [112131]$
- 计算 $A^{T} A$ ： $A^{T} A = [112131] 123111 = [1 \cdot 1 + 2 \cdot 2 + 3 \cdot 3 1 \cdot 1 + 1 \cdot 2 + 1 \cdot 3 1 \cdot 1 + 2 \cdot 1 + 3 \cdot 1 1 \cdot 1 + 1 \cdot 1 + 1 \cdot 1] = [14663]$
- 计算 $A^{T} b$ ： $A^{T} b = [112131] 235 = [1 \cdot 2 + 2 \cdot 3 + 3 \cdot 5 1 \cdot 2 + 1 \cdot 3 + 1 \cdot 5] = [2310]$
求解线性方程组：解 $A^{T} A x = A^{T} b$ ，即： $[14663] [a b] = [2310]$ 使用高斯消元法或矩阵求逆：
- 矩阵行列式： $14 \cdot 3 - 6 \cdot 6 = 42 - 36 = 6$
- 求逆矩阵： $(A^{T} A)^{- 1} = \frac{1}{6} [3 - 6 - 6 14] = [0.5 - 1 - 1 \frac{14}{6}] = [0.5 - 1 - 1 2.333...]$
- 计算 $x = (A^{T} A)^{- 1} A^{T} b$ ： $[a b] = [0.5 - 1 - 1 2.333...] [2310] = [0.5 \cdot 23 + (- 1) \cdot 10 (- 1) \cdot 23 + 2.333... \cdot 10] = [11.5 - 10 - 23 + 23.333...] = [1.5 0.333...]$ 因此， $a = 1.5$ ， $b = \frac{1}{3} \approx 0.333$ 。
结果：拟合直线为 $y = 1.5 x + 0.333$ 。可以通过计算残差 $∣∣ A x - b ∣ ∣^{2}$ 验证拟合效果。

3. 非线性最小二乘（Nonlinear Least Squares）

3.1 问题描述

拟合模型 $y = f (x, θ)$ ， $θ$ 为参数， $f$ 非线性
目标：最小化 $S (θ) = \sum_{i = 1}^{m} (y_{i} - f (x_{i}, θ))^{2}$

3.2 Gauss-Newton 方法

思想：用泰勒展开线性化，迭代求解
步骤：
1. 初始猜测 $θ_{0}$
2. 计算残差 $r_{i} = y_{i} - f (x_{i}, θ)$
3. 计算雅可比矩阵 $J_{ij} = \frac{\partial f ( x _{i} , θ )}{\partial θ _{j}}$
4. 解正规方程 $J^{T} J Δ θ = J^{T} r$
5. 更新 $θ \leftarrow θ + Δ θ$
6. 迭代至收敛

3.3 典型例题

例题2：用最小二乘法拟合二次曲线 $y = a x^{2} + b x + c$ ，数据点 $(1, 2), (2, 3), (3, 7), (4, 8)$

详细解答：

问题建模：虽然二次曲线是非线性函数，但对于参数 $a, b, c$ 来说，模型是线性的，因此可以用线性最小二乘法直接求解。每个数据点对应一个方程：
- 对于 $(1, 2)$ ： $a \cdot 1^{2} + b \cdot 1 + c = 2$
- 对于 $(2, 3)$ ： $a \cdot 2^{2} + b \cdot 2 + c = 3$
- 对于 $(3, 7)$ ： $a \cdot 3^{2} + b \cdot 3 + c = 7$
- 对于 $(4, 8)$ ： $a \cdot 4^{2} + b \cdot 4 + c = 8$
矩阵形式：将方程组写成矩阵形式 $A x = b$ ： $1491612341111 a b c = 2378$
正规方程：计算 $A^{T} A$ 和 $A^{T} b$ ：
- $A^{T} = 1114219311641$
- $A^{T} A = 1 + 16 + 81 + 256 1 + 8 + 27 + 64 1 + 4 + 9 + 16 1 + 8 + 27 + 64 1 + 4 + 9 + 16 1 + 2 + 3 + 4 1 + 4 + 9 + 16 1 + 2 + 3 + 4 1 + 1 + 1 + 1 = 35410030100301030104$
- $A^{T} b = 1 \cdot 2 + 4 \cdot 3 + 9 \cdot 7 + 16 \cdot 8 1 \cdot 2 + 2 \cdot 3 + 3 \cdot 7 + 4 \cdot 8 1 \cdot 2 + 1 \cdot 3 + 1 \cdot 7 + 1 \cdot 8 = 2056120$
求解线性方程组：解 $A^{T} A x = A^{T} b$ ，即： $35410030100301030104 a b c = 2056120$ 使用数值方法或矩阵求逆，得到：
- $a \approx 0.000$ （几乎为0）
- $b \approx 2.200$
- $c \approx - 0.500$
结果：拟合曲线为 $y = 2.2 x - 0.5$ 。这表明数据点更适合用直线而非二次曲线进行拟合。可以通过计算残差验证拟合效果。

例题3：用非线性最小二乘法拟合指数模型 $y = a e^{b x}$ ，数据点 $(0, 1), (1, 2), (2, 5)$

详细解答：

问题建模：目标是找到参数 $a$ 和 $b$ ，使得 $y = a e^{b x}$ 拟合数据点。模型为非线性，无法直接用线性最小二乘法求解，因此使用Gauss-Newton迭代方法。
目标函数：最小化残差平方和： $S (a, b) = i = 1 \sum 3 (y_{i} - a e^{b x_{i}})^{2}$
初始猜测：选择初始值 $a = 1, b = 0.5$ （可以根据数据趋势粗略估计）。
迭代步骤：
- 计算残差：对于当前参数 $a, b$ ，计算每个数据点的残差 $r_{i} = y_{i} - a e^{b x_{i}}$ 。
  - $x = 0$ ： $r_{1} = 1 - 1 \cdot e^{0.5 \cdot 0} = 1 - 1 = 0$
  - $x = 1$ ： $r_{2} = 2 - 1 \cdot e^{0.5 \cdot 1} \approx 2 - 1.6487 = 0.3513$
  - $x = 2$ ： $r_{3} = 5 - 1 \cdot e^{0.5 \cdot 2} \approx 5 - 2.7183 = 2.2817$
- 计算雅可比矩阵： $J_{ij} = \frac{\partial f ( x _{i} , θ )}{\partial θ _{j}}$ ，其中 $f (x, a, b) = a e^{b x}$ ， $θ_{1} = a$ ， $θ_{2} = b$ 。
  - $\frac{\partial f}{\partial a} = e^{b x}$
  - $\frac{\partial f}{\partial b} = a x e^{b x}$
  - 对于 $x = 0$ ： $J_{11} = e^{0.5 \cdot 0} = 1$ ， $J_{12} = 0 \cdot 1 \cdot e^{0.5 \cdot 0} = 0$
  - 对于 $x = 1$ ： $J_{21} = e^{0.5 \cdot 1} \approx 1.6487$ ， $J_{22} = 1 \cdot 1 \cdot e^{0.5 \cdot 1} \approx 1.6487$
  - 对于 $x = 2$ ： $J_{31} = e^{0.5 \cdot 2} \approx 2.7183$ ， $J_{32} = 2 \cdot 1 \cdot e^{0.5 \cdot 2} \approx 5.4366$
  - 雅可比矩阵： $J = 1 1.6487 2.7183 0 1.6487 5.4366$
- 解正规方程：计算 $J^{T} J$ 和 $J^{T} r$ ，解 $J^{T} J Δ θ = J^{T} r$ ，得到参数更新量 $Δ a, Δ b$ 。
- 更新参数： $a \leftarrow a + Δ a$ ， $b \leftarrow b + Δ b$ 。
- 重复迭代：重复上述步骤，直到参数收敛或残差足够小。
结果：经过多次迭代，参数收敛到 $a \approx 0.884$ ， $b \approx 0.864$ ，拟合模型为 $y \approx 0.884 e^{0.864 x}$ ，接近真实指数增长趋势。

4. GPS定位中的最小二乘

原理：通过测量与多颗卫星的距离，建立超定方程组，利用最小二乘法估算接收机位置。
模型： $∣∣ x - s_{i} ∣∣ = d_{i}$ ， $x$ 为未知位置， $s_{i}$ 为卫星位置， $d_{i}$ 为测距
步骤：线性化后用Gauss-Newton迭代

4.1 典型例题

例题4：假设接收机与三颗卫星的距离测量如下，估算接收机位置 $(x, y)$ ：

卫星1：位置 $(0, 0)$ ，距离 $d_{1} = 5$
卫星2：位置 $(8, 0)$ ，距离 $d_{2} = 5$
卫星3：位置 $(4, 6)$ ，距离 $d_{3} = 5$

详细解答：

问题建模：接收机位置为 $(x, y)$ ，与各卫星的距离方程为：
- $(x - 0)^{2} + (y - 0)^{2} = 5 ⟹ x^{2} + y^{2} = 25$
- $(x - 8)^{2} + (y - 0)^{2} = 5 ⟹ (x - 8)^{2} + y^{2} = 25$
- $(x - 4)^{2} + (y - 6)^{2} = 5 ⟹ (x - 4)^{2} + (y - 6)^{2} = 25$
非线性方程组：上述方程是非线性的，直接求解较复杂，因此使用Gauss-Newton方法进行迭代求解。
初始猜测：根据几何图形，接收机可能位于三颗卫星的中心附近，初始猜测为 $(x, y) = (4, 2)$ 。
线性化与迭代：
- 定义残差函数：
  - $r_{1} = x^{2} + y^{2} - 5$
  - $r_{2} = (x - 8)^{2} + y^{2} - 5$
  - $r_{3} = (x - 4)^{2} + (y - 6)^{2} - 5$
- 计算雅可比矩阵（偏导数）：
  - $\frac{\partial r _{1}}{\partial x} = \frac{x}{x ^{2} + y ^{2}}$ ， $\frac{\partial r _{1}}{\partial y} = \frac{y}{x ^{2} + y ^{2}}$
  - $\frac{\partial r _{2}}{\partial x} = \frac{x - 8}{( x - 8 ) ^{2} + y ^{2}}$ ， $\frac{\partial r _{2}}{\partial y} = \frac{y}{( x - 8 ) ^{2} + y ^{2}}$
  - $\frac{\partial r _{3}}{\partial x} = \frac{x - 4}{( x - 4 ) ^{2} + ( y - 6 ) ^{2}}$ ， $\frac{\partial r _{3}}{\partial y} = \frac{y - 6}{( x - 4 ) ^{2} + ( y - 6 ) ^{2}}$
- 在初始点 $(4, 2)$ 处计算残差和雅可比矩阵，解正规方程更新参数。
- 经过多次迭代，位置收敛到 $(x, y) \approx (4, 1.5)$ 。
结果：接收机位置为 $(4, 1.5)$ ，可以通过代入原方程验证各卫星到该位置的距离是否接近5。

5. 常见考点与易错点

正规方程推导与矩阵乘法
QR分解步骤与正交化方法（Gram-Schmidt）
非线性最小二乘的迭代过程
残差（residual）、条件数（condition number）理解
英文术语拼写与公式记忆

6. 英文关键词与表达

least squares, overdetermined system, normal equation, QR decomposition, orthogonalization, Gram-Schmidt, residual, regression, nonlinear least squares, Gauss-Newton, GPS positioning, condition number

如需更详细例题推导或某一算法的代码实现，可随时补充！

CE7453 Numerical Algorithms 期末高分超详细攻略（第六章：Eigenanalysis）

本章内容：特征值与特征向量、幂迭代、QR算法、PageRank
适用对象：零基础/考前冲刺/快速查漏补缺
内容特色：详细原理、公式推导、算法流程、例题全解、英文关键词

1. 基本概念与背景

特征值与特征向量（Eigenvalue & Eigenvector）：对于 $n \times n$ 矩阵 $A$ ，如果存在非零向量 $x$ 和标量 $λ$ 使 $A x = λ x$ ，则 $λ$ 为特征值， $x$ 为对应特征向量。
实际应用：PCA主成分分析、PageRank、振动分析、图像处理、机器学习等。
数学意义：特征值表示矩阵在某些方向上的缩放因子，特征向量表示这些方向。特征分解可以帮助理解矩阵的性质和行为。

特征值分解（Eigenvalue Decomposition）：如果矩阵 $A$ 可对角化，则可以表示为 $A = Q Λ Q^{- 1}$ ，其中 $Q$ 是特征向量组成的矩阵， $Λ$ 是特征值组成对角矩阵。这种分解在求矩阵幂、解微分方程等场景中非常有用。

2. 幂迭代法（Power Iteration）

2.1 原理

用于求最大模特征值（dominant eigenvalue）及其对应的特征向量。
基本思想是通过反复用矩阵 $A$ 作用于初始向量，逐步放大主特征分量（与最大特征值对应的分量），最终收敛到主特征向量。
数学依据：假设矩阵 $A$ 有特征值 $λ_{1}, λ_{2}, ..., λ_{n}$ ，且 $∣ λ_{1} ∣ > ∣ λ_{2} ∣ \geq ... \geq ∣ λ_{n} ∣$ ，初始向量 $x_{0}$ 可以表示为特征向量的线性组合 $x_{0} = c_{1} v_{1} + c_{2} v_{2} + ... + c_{n} v_{n}$ ，则 $A^{k} x_{0} \approx c_{1} λ_{1}^{k} v_{1}$ （当 $k$ 很大时），即收敛到主特征向量方向。

2.2 算法流程

选择一个非零初始向量 $x_{0}$ （通常随机选择或全为1）。
归一化 $x_{0}$ ，以避免数值过大或过小。
迭代：计算 $x_{k + 1} = A x_{k}$ 。
归一化 $x_{k + 1}$ ，得到单位向量。
检查收敛性：若 $∣∣ x_{k + 1} - x_{k} ∣∣ < ϵ$ （或达到最大迭代次数），则停止迭代。
计算特征值近似： $λ \approx \frac{x _{k}^{T} A x _{k}}{x _{k}^{T} x _{k}}$ （Rayleigh quotient，瑞利商）。

注意事项：

初始向量不能与主特征向量正交，否则无法收敛到最大特征值。
收敛速度取决于特征值之间的比值 $∣ λ_{2} / λ_{1} ∣$ ，比值越小，收敛越快。

2.3 例题

例题1：给定矩阵 $A = [2112]$ ，初始向量 $x_{0} = [10]$ ，用幂迭代法求最大特征值及其特征向量。

详细解答：

初始化： $x_{0} = [10]$ ，归一化后仍为 $[10]$ （因为 $∣∣ x_{0} ∣ ∣_{2} = 1$ ）。
第一次迭代：
- 计算 $x_{1} = A x_{0} = [2112] [10] = [21]$ 。
- 归一化： $∣∣ x_{1} ∣ ∣_{2} = 2^{2} + 1^{2} = 5 \approx 2.236$ ，所以 $x_{1} = [2/2.236 1/2.236] \approx [0.894 0.447]$ 。
第二次迭代：
- 计算 $x_{2} = A x_{1} = [2112] [0.894 0.447] = [2 * 0.894 + 1 * 0.447 1 * 0.894 + 2 * 0.447] = [2.235 1.788]$ 。
- 归一化： $∣∣ x_{2} ∣ ∣_{2} = 2.23 5^{2} + 1.78 8^{2} \approx 2.86$ ，所以 $x_{2} \approx [0.781 0.625]$ 。
第三次迭代：
- 计算 $x_{3} = A x_{2} \approx [2 * 0.781 + 1 * 0.625 1 * 0.781 + 2 * 0.625] = [2.187 2.031]$ 。
- 归一化： $∣∣ x_{3} ∣ ∣_{2} \approx 2.98$ ，所以 $x_{3} \approx [0.734 0.681]$ 。
第四次迭代：
- 计算 $x_{4} = A x_{3} \approx [2 * 0.734 + 1 * 0.681 1 * 0.734 + 2 * 0.681] = [2.149 2.096]$ 。
- 归一化： $∣∣ x_{4} ∣ ∣_{2} \approx 3.00$ ，所以 $x_{4} \approx [0.716 0.698]$ 。
第五次迭代：
- 计算 $x_{5} = A x_{4} \approx [2 * 0.716 + 1 * 0.698 1 * 0.716 + 2 * 0.698] = [2.130 2.112]$ 。
- 归一化： $∣∣ x_{5} ∣ ∣_{2} \approx 3.00$ ，所以 $x_{5} \approx [0.710 0.704]$ 。
收敛检查：继续迭代，向量逐渐接近 $[0.707 0.707]$ （即 $[11]$ 归一化后），说明收敛到主特征向量。
计算特征值：使用瑞利商， $λ \approx \frac{x _{5}^{T} A x _{5}}{x _{5}^{T} x _{5}}$ 。
- $A x_{5} \approx [2 * 0.710 + 1 * 0.704 1 * 0.710 + 2 * 0.704] = [2.124 2.118]$ 。
- $x_{5}^{T} A x_{5} \approx 0.710 * 2.124 + 0.704 * 2.118 \approx 3.000$ 。
- $x_{5}^{T} x_{5} = 0.71 0^{2} + 0.70 4^{2} \approx 1$ 。
- 因此 $λ \approx 3$ 。

答案：最大特征值约为 $λ = 3$ ，对应特征向量为 $[11]$ （未归一化形式）。

验证：直接计算特征值， $A$ 的特征方程为 $det (A - λ I) = (2 - λ)^{2} - 1 = 0$ ，解得 $λ = 3$ 或 $λ = 1$ ，与幂迭代结果一致。

3. QR算法（QR Algorithm）

3.1 原理

QR算法是一种求解矩阵所有特征值的迭代方法，特别适合对称矩阵。
基本思想：通过反复对矩阵进行QR分解（将矩阵分解为正交矩阵 $Q$ 和上三角矩阵 $R$ 的乘积），并重新组合为 $A_{k + 1} = R_{k} Q_{k}$ ，最终使矩阵收敛到对角形式（或近似对角形式），对角线元素即为特征值。
数学依据：QR算法本质上是幂法的扩展，同时对所有特征值进行迭代，收敛后矩阵变为上三角或对角矩阵。

3.2 算法流程

初始化：令 $A_{0} = A$ 。
对当前矩阵 $A_{k}$ 进行QR分解： $A_{k} = Q_{k} R_{k}$ （其中 $Q_{k}$ 是正交矩阵， $R_{k}$ 是上三角矩阵）。
计算新的矩阵： $A_{k + 1} = R_{k} Q_{k}$ 。
重复步骤2和3，直到 $A_{k}$ 近似为对角矩阵（或达到最大迭代次数），对角线上的元素即为特征值的近似值。

注意事项：

对于对称矩阵，QR算法会收敛到对角矩阵；对于非对称矩阵，收敛到上三角矩阵（Schur形式）。
实际应用中常结合Hessenberg化（将矩阵化为上Hessenberg形式）和位移策略（shift）来加速收敛。

3.3 例题

例题2：使用QR算法计算矩阵 $A = [5445]$ 的特征值。

详细解答：

初始化： $A_{0} = [5445]$ 。
第一次迭代：
- 对 $A_{0}$ 进行QR分解：
  - 使用Gram-Schmidt正交化或其他方法，得到 $Q_{0} = [0.707 0.707 - 0.707 0.707]$ ， $R_{0} = [7.071 0 6.364 1.414]$ （近似值）。
- 计算 $A_{1} = R_{0} Q_{0} = [7.071 0 6.364 1.414] [0.707 0.707 - 0.707 0.707] \approx [9.5 1.0 - 0.5 0.5]$ 。
第二次迭代：
- 对 $A_{1}$ 进行QR分解，得到新的 $Q_{1}$ 和 $R_{1}$ 。
- 计算 $A_{2} = R_{1} Q_{1}$ ，结果更接近对角矩阵。
继续迭代：经过多次迭代，矩阵逐渐收敛到 $[9001]$ ，即特征值为 $9$ 和 $1$ 。

答案：特征值为 $λ_{1} = 9$ ， $λ_{2} = 1$ 。

验证：直接计算特征方程 $det (A - λ I) = (5 - λ)^{2} - 16 = 0$ ，解得 $λ = 9$ 或 $λ = 1$ ，与QR算法结果一致。

4. PageRank 算法（Google 搜索排序）

4.1 背景

PageRank 是Google搜索引擎的核心算法之一，用于评估网页的重要性。
基本思想：网页的重要性取决于链接到该网页的其他网页数量和重要性，本质上是求网页链接矩阵的主特征向量。
数学模型：设 $P$ 为转移概率矩阵（网页之间的链接概率），PageRank向量 $r$ 满足 $P r = r$ ，即 $r$ 是矩阵 $P$ 对应特征值 $λ = 1$ 的特征向量。

4.2 算法流程

构造转移矩阵 $P$ ： $P_{ij}$ 表示从网页 $j$ 链接到网页 $i$ 的概率（若网页 $j$ 有 $k$ 个外链，则每个外链概率为 $1/ k$ ）。
选择初始向量 $r_{0}$ （通常为均匀分布，即每个网页初始重要性相等）。
迭代： $r_{k + 1} = P r_{k}$ 。
归一化 $r_{k + 1}$ ，确保向量元素之和为1。
收敛后， $r$ 的各个分量即为各网页的PageRank值（排名分数）。

阻尼因子（Damping Factor）：

实际中引入阻尼因子 $d$ （通常为0.85），以模拟用户随机跳转行为，修正公式为 $r = (1 - d) / n + d P r$ ，其中 $n$ 为网页总数。
这保证了矩阵的收敛性，避免某些网页没有外链导致的问题。

4.3 例题

例题3：假设有3个网页，链接关系如下：网页1链接到2和3，网页2链接到3，网页3链接到1。计算各网页的PageRank值（忽略阻尼因子）。

详细解答：

构造转移矩阵 $P$ ：
- 网页1有2个外链（到2和3），所以 $P_{21} = 0.5$ ， $P_{31} = 0.5$ ， $P_{11} = 0$ 。
- 网页2有1个外链（到3），所以 $P_{32} = 1$ ， $P_{12} = 0$ ， $P_{22} = 0$ 。
- 网页3有1个外链（到1），所以 $P_{13} = 1$ ， $P_{23} = 0$ ， $P_{33} = 0$ 。
- 转移矩阵 $P = 0 0.5 0.5 001100$ 。
初始化： $r_{0} = 1/3 1/3 1/3$ 。
第一次迭代：
- $r_{1} = P r_{0} = 0 0.5 0.5 001100 1/3 1/3 1/3 = 1/3 1/6 1/2$ 。
- 归一化（可选，此处不归一化以观察收敛）。
第二次迭代：
- $r_{2} = P r_{1} = 0 0.5 0.5 001100 1/3 1/6 1/2 = 1/2 1/6 1/3$ 。
继续迭代：经过多次迭代， $r$ 收敛到 $0.4 0.2 0.4$ （近似值）。

答案：网页1的PageRank值为 $0.4$ ，网页2的值为 $0.2$ ，网页3的值为 $0.4$ ，因此网页1和网页3最重要且同等重要。

5. 修正例题：幂迭代法

例题4（修正新增例题3）：计算矩阵 $A = [4113]$ 的最大特征值及其特征向量。

详细解答：

初始化： $x_{0} = [10]$ ，归一化后仍为 $[10]$ 。
第一次迭代：
- $x_{1} = A x_{0} = [4113] [10] = [41]$ 。
- 归一化： $∣∣ x_{1} ∣ ∣_{2} = 16 + 1 = 17 \approx 4.123$ ， $x_{1} \approx [0.970 0.243]$ 。
第二次迭代：
- $x_{2} = A x_{1} \approx [4 * 0.970 + 1 * 0.243 1 * 0.970 + 3 * 0.243] = [4.123 1.699]$ 。
- 归一化： $∣∣ x_{2} ∣ ∣_{2} \approx 4.456$ ， $x_{2} \approx [0.925 0.381]$ 。
第三次迭代：
- $x_{3} = A x_{2} \approx [4 * 0.925 + 1 * 0.381 1 * 0.925 + 3 * 0.381] = [4.081 2.068]$ 。
- 归一化： $∣∣ x_{3} ∣ ∣_{2} \approx 4.564$ ， $x_{3} \approx [0.894 0.453]$ 。
继续迭代：向量逐渐收敛到 $[0.850 0.527]$ ，即未归一化形式为 $[1.615 1]$ 或近似 $[1.618 1]$ 。
计算特征值： $λ \approx \frac{x _{k}^{T} A x _{k}}{x _{k}^{T} x _{k}}$ ，取 $x_{k} \approx [0.850 0.527]$ ：
- $A x_{k} \approx [4 * 0.850 + 1 * 0.527 1 * 0.850 + 3 * 0.527] = [3.927 2.431]$ 。
- $x_{k}^{T} A x_{k} \approx 0.850 * 3.927 + 0.527 * 2.431 \approx 4.618$ 。
- $x_{k}^{T} x_{k} = 0.85 0^{2} + 0.52 7^{2} \approx 1$ 。
- 因此 $λ \approx 4.618$ 。

答案：最大特征值约为 $λ \approx 4.618$ ，对应特征向量为 $[1.618 1]$ （近似值）。

验证：特征方程 $det (A - λ I) = (4 - λ) (3 - λ) - 1 = λ^{2} - 7 λ + 11 = 0$ ，解得 $λ = \frac{7 \pm 5}{2}$ ，即 $λ \approx 4.618$ 或 $λ \approx 2.382$ ，与幂迭代结果一致。

6. 常见考点与易错点

幂迭代法：初始向量不能与主特征向量正交，否则无法收敛到最大特征值；收敛速度受特征值比值影响。
QR算法：只适合小规模矩阵或对称矩阵，实际应用中需结合位移策略；理解QR分解的正交性和收敛原理。
特征值/向量定义与实际意义：特征值分解在矩阵分析和应用中的重要性。
PageRank：转移矩阵的构造方法，阻尼因子的作用，归一化的必要性。
英文术语拼写与公式记忆：确保考试中能准确写出专业术语和公式。

7. 英文关键词与表达

eigenvalue, eigenvector, power iteration, Rayleigh quotient, QR algorithm, diagonalization, convergence, PageRank, principal component analysis (PCA), dominant eigenvalue, damping factor, transition matrix

如需更详细例题推导或某一算法的代码实现，可随时补充！

CE7453 Numerical Algorithms 期末高分超详细攻略（第七章：Fourier Transform）

本章内容：离散傅里叶变换（DFT）、快速傅里叶变换（FFT）、三角插值、JPEG应用
适用对象：零基础/考前冲刺/快速查漏补缺
内容特色：详细原理、公式推导、算法流程、例题全解、英文关键词

1. 基本概念与背景

傅里叶变换（Fourier Transform）：将信号从时域变换到频域，分析信号的频率成分。
离散傅里叶变换（DFT）：对有限长离散信号进行傅里叶变换，常用于数字信号处理、图像压缩等。
快速傅里叶变换（FFT）：高效计算DFT的算法，大幅降低运算量。
实际应用：音频分析、图像处理（如JPEG压缩）、通信、工程仿真等。

2. 离散傅里叶变换（DFT）

2.1 公式与推导

DFT公式：
$X_{k} = n = 0 \sum N - 1 x_{n} e^{- 2 πikn / N}, k = 0, 1, ..., N - 1$ 其中， $X_{k}$ 是频域信号， $x_{n}$ 是时域信号， $N$ 是信号长度， $i$ 是虚数单位。
IDFT公式（逆变换）：
$x_{n} = \frac{1}{N} k = 0 \sum N - 1 X_{k} e^{2 πikn / N}$ IDFT 将频域信号转换回时域信号。
性质：
- 可逆性：DFT 和 IDFT 互为逆运算，可以无损地从时域到频域再回到时域。
- 能量守恒：根据 Parseval 定理，信号在时域和频域的能量相等。
- 周期性：DFT 结果具有周期性， $X_{k + N} = X_{k}$ 。

2.2 计算复杂度

直接计算 DFT 的复杂度为 $O (N^{2})$ ，因为对于每个 $k$ ，需要对所有 $n$ 进行求和运算。
当信号长度 $N$ 很大时，计算效率极低，因此需要更高效的算法（如 FFT）。

3. 快速傅里叶变换（FFT）

3.1 原理

分治思想：将长度为 $N$ 的 DFT 分解为两个长度为 $N /2$ 的 DFT，递归进行分解，直到长度为 1。
要求： $N$ 必须为 2 的幂，以便递归分解。
FFT 通过减少重复计算，大幅降低复杂度。

3.2 算法流程（Cooley-Tukey FFT）

输入序列 $x_{0}, x_{1}, ..., x_{N - 1}$ 。
若 $N = 1$ ，直接返回 $x_{0}$ 。
将序列分为偶数项 $x_{2 m}$ 和奇数项 $x_{2 m + 1}$ 。
分别递归计算偶数和奇数部分的 DFT，得到 $E_{k}$ 和 $O_{k}$ 。
合并结果：
$X_{k} = E_{k} + W_{N}^{k} O_{k} X_{k + N /2} = E_{k} - W_{N}^{k} O_{k}$ 其中 $E_{k}$ 、 $O_{k}$ 分别为偶/奇部分的 DFT 结果， $W_{N}^{k} = e^{- 2 πik / N}$ 是旋转因子。
复杂度从 $O (N^{2})$ 降为 $O (N lo g N)$ 。

3.3 例题

例题1：计算 $x = [1, 0, - 1, 0]$ 的 DFT。

详细解答：

给定 $N = 4$ ，使用 DFT 公式 $X_{k} = \sum_{n = 0}^{3} x_{n} e^{- 2 πikn /4}$ 。
计算 $X_{0}$ ： $X_{0} = x_{0} e^{0} + x_{1} e^{0} + x_{2} e^{0} + x_{3} e^{0} = 1 + 0 + (- 1) + 0 = 0$
计算 $X_{1}$ ： $X_{1} = x_{0} e^{0} + x_{1} e^{- 2 πi \cdot 1 \cdot 1/4} + x_{2} e^{- 2 πi \cdot 2 \cdot 1/4} + x_{3} e^{- 2 πi \cdot 3 \cdot 1/4}$ 其中 $e^{- 2 πi /4} = e^{- πi /2} = - i$ ， $e^{- πi} = - 1$ ， $e^{- 3 πi /2} = i$ ，所以： $X_{1} = 1 \cdot 1 + 0 \cdot (- i) + (- 1) \cdot (- 1) + 0 \cdot i = 1 + 0 + 1 + 0 = 2$
计算 $X_{2}$ ： $X_{2} = x_{0} e^{0} + x_{1} e^{- 2 πi \cdot 1 \cdot 2/4} + x_{2} e^{- 2 πi \cdot 2 \cdot 2/4} + x_{3} e^{- 2 πi \cdot 3 \cdot 2/4}$ 其中 $e^{- πi} = - 1$ ， $e^{- 2 πi} = 1$ ， $e^{- 3 πi} = - 1$ ，所以： $X_{2} = 1 \cdot 1 + 0 \cdot (- 1) + (- 1) \cdot 1 + 0 \cdot (- 1) = 1 + 0 + (- 1) + 0 = 0$
计算 $X_{3}$ ： $X_{3} = x_{0} e^{0} + x_{1} e^{- 2 πi \cdot 1 \cdot 3/4} + x_{2} e^{- 2 πi \cdot 2 \cdot 3/4} + x_{3} e^{- 2 πi \cdot 3 \cdot 3/4}$ 其中 $e^{- 3 πi /2} = i$ ， $e^{- 3 πi} = - 1$ ， $e^{- 9 πi /2} = e^{- 4 πi - πi /2} = e^{- πi /2} = - i$ ，所以： $X_{3} = 1 \cdot 1 + 0 \cdot i + (- 1) \cdot (- 1) + 0 \cdot (- i) = 1 + 0 + 1 + 0 = 2$
最终结果： $X = [0, 2, 0, 2]$ 。

例题2：使用 FFT 计算 $x = [1, 2, 3, 4, 5, 6, 7, 8]$ 的 DFT。

解答步骤：

$N = 8$ ，满足 2 的幂要求。
第一步分解：分为偶数项 $[1, 3, 5, 7]$ 和奇数项 $[2, 4, 6, 8]$ 。
继续递归分解，直到长度为 1。
合并计算（由于篇幅限制，这里省略详细计算步骤，直接给出结果）。
结果： $X = [36, - 4 + 9.656 i, - 4 + 4 i, - 4 + 1.656 i, - 4, - 4 - 1.656 i, - 4 - 4 i, - 4 - 9.656 i]$ 。

4. 三角插值与JPEG应用

4.1 三角插值（Trigonometric Interpolation）

定义：用三角函数（正弦、余弦）拟合周期数据，适用于周期性信号的分析和重构。
与 DFT 的关系：DFT 本质上就是三角插值系数的计算， $X_{k}$ 表示信号在不同频率上的分量。
应用：信号重构、数据平滑、周期性预测。

例题3：计算信号 $x = [1, 1, 1, 1]$ 的 DFT，并解释其频域意义。

详细解答：

使用 DFT 公式： $X_{k} = \sum_{n = 0}^{3} x_{n} e^{- 2 πikn /4}$ 。
计算 $X_{0}$ ： $X_{0} = 1 + 1 + 1 + 1 = 4$
计算 $X_{1}$ ： $X_{1} = 1 \cdot 1 + 1 \cdot e^{- πi /2} + 1 \cdot e^{- πi} + 1 \cdot e^{- 3 πi /2} = 1 + (- i) + (- 1) + i = 0$
计算 $X_{2}$ ： $X_{2} = 1 \cdot 1 + 1 \cdot e^{- πi} + 1 \cdot e^{- 2 πi} + 1 \cdot e^{- 3 πi} = 1 + (- 1) + 1 + (- 1) = 0$
计算 $X_{3}$ ： $X_{3} = 1 \cdot 1 + 1 \cdot e^{- 3 πi /2} + 1 \cdot e^{- 3 πi} + 1 \cdot e^{- 9 πi /2} = 1 + i + (- 1) + (- i) = 0$
结果： $X = [4, 0, 0, 0]$ 。
频域意义：信号的能量全部集中在直流分量（ $X_{0}$ ），表示信号没有周期性变化，是一个常数值信号。

4.2 JPEG 压缩中的 DCT

离散余弦变换（DCT）：JPEG 采用 DCT，本质与 DFT 类似，但只使用实数部分（余弦函数），避免复数运算。
DCT 公式（一维）： $C_{k} = n = 0 \sum N - 1 x_{n} cos [\frac{π}{N} (n + \frac{1}{2}) k]$
优势：
- DCT 能将图像能量集中在低频分量，便于去除高频噪声（人眼对高频细节不敏感）。
- 实数运算比 DFT 的复数运算更高效。
JPEG 压缩流程：
1. 将图像分成 8x8 像素块。
2. 对每个块应用二维 DCT，得到频域系数。
3. 量化：对高频系数进行压缩（丢弃小值或近似）。
4. 编码：使用霍夫曼编码等方法进一步压缩数据。
5. 解压时逆向操作：逆量化、逆 DCT 重构图像。
DCT 与 DFT 的区别：
- DCT 仅用余弦函数，输出为实数；DFT 使用复指数，输出为复数。
- DCT 边界条件更适合图像处理（镜像对称），减少边界不连续性。

例题4：解释为什么 JPEG 压缩会丢失细节。

解答：

JPEG 压缩通过量化步骤丢弃高频分量的小值系数，而高频分量对应图像的细节（如边缘、纹理）。
人眼对高频细节不敏感，因此丢弃这些分量可以在不明显影响视觉质量的情况下大幅减少数据量。
这种有损压缩导致细节丢失，尤其在高压缩比下，图像可能出现块状伪影（blocking artifacts）。

5. 常见考点与易错点

DFT/IDFT 公式推导与正负号：注意指数中的正负号，DFT 为负号，IDFT 为正号。
FFT 分治递归与合并步骤：理解如何分解序列和合并结果，旋转因子的作用。
频域与时域的物理意义：时域表示信号随时间变化，频域表示信号的频率成分和能量分布。
DCT 与 DFT 的区别：DCT 仅用实数，适合图像压缩；DFT 处理复数，适合信号分析。
易错点：
- 计算 DFT 时漏掉旋转因子的复数性质。
- FFT 合并时未正确处理 $N /2$ 偏移。
- 混淆 DCT 和 DFT 的应用场景。
英文术语拼写与公式记忆：熟悉常见术语拼写，考试中可能要求手写公式。

6. 英文关键词与表达

discrete Fourier transform (DFT), fast Fourier transform (FFT), frequency domain, time domain, trigonometric interpolation, discrete cosine transform (DCT), JPEG compression, signal processing, spectrum, complex exponential, quantization, blocking artifacts

如需更详细例题推导或某一算法的代码实现，可随时补充！

CE7453 Numerical Algorithms 期末高分超详细攻略（第八章：3D Data Registration）

本章内容：3D数据配准、ICP算法、最优旋转/平移、SVD分解、实际应用
适用对象：零基础/考前冲刺/快速查漏补缺
内容特色：详细原理、公式推导、算法流程、例题全解、英文关键词

1. 基本概念与背景

3D数据配准（3D Data Registration）：将多个3D点云或模型对齐到同一坐标系，常用于3D重建、医学成像、机器人导航等。
典型问题：给定两组点云，求最优刚性变换（旋转+平移）使它们尽量重合。

2. ICP算法（Iterative Closest Point）

2.1 原理

迭代式地将源点云与目标点云配准
每次迭代包括“最近点匹配”和“最优变换估计”两步

2.2 算法详细流程

初始化：设源点集 $A$ ，目标点集 $B$ ，初始变换 $T_{0}$
最近点匹配：对 $A$ 中每个点 $a_{i}$ ，在 $B$ 中找最近点 $b_{j}$
估计最优变换：求刚性变换 $T$ （旋转 $R$ +平移 $t$ ），最小化配对点的均方误差 $R, t min i \sum ∥ R a_{i} + t - b_{i} ∥^{2}$
应用变换：用 $T$ 更新 $A$
收敛判定：若误差变化小于阈值，停止；否则回到第2步

2.3 最优旋转/平移的SVD解法

设 $A$ 和 $B$ 已配对，先去中心化
计算协方差矩阵 $H = \sum_{i} (a_{i} - \overset{a}{ˉ}) (b_{i} - \overset{ˉ}{b})^{T}$
对 $H$ 做SVD分解 $H = U Σ V^{T}$
最优旋转 $R = V U^{T}$ ，最优平移 $t = \overset{ˉ}{b} - R \overset{a}{ˉ}$

3. 典型例题

例题1：已知 $A = {(0, 0), (1, 0)}$ ， $B = {(1, 1), (2, 1)}$ ，求将 $A$ 配准到 $B$ 的最优刚性变换

解答：

计算质心：
- $A$ 的质心： $\overset{a}{ˉ} = (\frac{0 + 1}{2}, \frac{0 + 0}{2}) = (0.5, 0)$
- $B$ 的质心： $\overset{ˉ}{b} = (\frac{1 + 2}{2}, \frac{1 + 1}{2}) = (1.5, 1)$
去中心化：
- $A$ 去中心化后： ${(0 - 0.5, 0 - 0), (1 - 0.5, 0 - 0)} = {(- 0.5, 0), (0.5, 0)}$
- $B$ 去中心化后： ${(1 - 1.5, 1 - 1), (2 - 1.5, 1 - 1)} = {(- 0.5, 0), (0.5, 0)}$
计算协方差矩阵 $H$ ：
- $H = \sum_{i} (a_{i} - \overset{a}{ˉ}) (b_{i} - \overset{ˉ}{b})^{T}$
- $H = (- 0.5) \cdot (- 0.5) + (0.5) \cdot (0.5) = 0.25 + 0.25 = 0.5$ （对每个维度分别计算）
- 完整矩阵 $H = [0.5 0 00]$
SVD分解：
- $H = U Σ V^{T}$ ，得到 $U = I$ ， $V = I$ （因为 $H$ 是对角矩阵）
- 旋转矩阵 $R = V U^{T} = I$ （单位矩阵）
计算平移向量：
- $t = \overset{ˉ}{b} - R \overset{a}{ˉ} = (1.5, 1) - (0.5, 0) = (1, 1)$
结果：
- 最优旋转矩阵 $R$ 为单位矩阵（无旋转）
- 最优平移向量 $t = (1, 1)$

验证：将 $A$ 的点应用变换后， $(0, 0) \to (1, 1)$ ， $(1, 0) \to (2, 1)$ ，与 $B$ 完全重合。

例题2：给定点集 $A = [(0, 0, 0), (1, 0, 0)]$ 和 $B = [(0, 1, 0), (1, 1, 0)]$ ，求将 $A$ 配准到 $B$ 的最优刚性变换

解答步骤：

计算质心：
- $A$ 的质心： $\overset{a}{ˉ} = (\frac{0 + 1}{2}, \frac{0 + 0}{2}, \frac{0 + 0}{2}) = (0.5, 0, 0)$
- $B$ 的质心： $\overset{ˉ}{b} = (\frac{0 + 1}{2}, \frac{1 + 1}{2}, \frac{0 + 0}{2}) = (0.5, 1, 0)$
去中心化：
- $A$ 去中心化后： ${(0 - 0.5, 0 - 0, 0 - 0), (1 - 0.5, 0 - 0, 0 - 0)} = {(- 0.5, 0, 0), (0.5, 0, 0)}$
- $B$ 去中心化后： ${(0 - 0.5, 1 - 1, 0 - 0), (1 - 0.5, 1 - 1, 0 - 0)} = {(- 0.5, 0, 0), (0.5, 0, 0)}$
计算协方差矩阵 $H$ ：
- $H = \sum_{i} (a_{i} - \overset{a}{ˉ}) (b_{i} - \overset{ˉ}{b})^{T}$
- 对于每个点对，计算外积并求和
- $H = 0.5 00 000000$
SVD分解：
- $H = U Σ V^{T}$ ，得到 $U = I$ ， $V = I$
- 旋转矩阵 $R = V U^{T} = I$ （单位矩阵）
计算平移向量：
- $t = \overset{ˉ}{b} - R \overset{a}{ˉ} = (0.5, 1, 0) - (0.5, 0, 0) = (0, 1, 0)$
结果：
- 最优旋转矩阵 $R$ 为单位矩阵（无旋转）
- 最优平移向量 $t = (0, 1, 0)$

验证：将 $A$ 的点应用变换后， $(0, 0, 0) \to (0, 1, 0)$ ， $(1, 0, 0) \to (1, 1, 0)$ ，与 $B$ 完全重合。

例题3（进阶）：给定点集 $A = [(0, 0, 0), (1, 0, 0), (0, 1, 0)]$ 和 $B = [(0, 0, 0), (0, 1, 0), (1, 0, 0)]$ ，求将 $A$ 配准到 $B$ 的最优刚性变换（涉及旋转）。

解答步骤：

计算质心：
- $A$ 的质心： $\overset{a}{ˉ} = (\frac{0 + 1 + 0}{3}, \frac{0 + 0 + 1}{3}, \frac{0 + 0 + 0}{3}) = (\frac{1}{3}, \frac{1}{3}, 0)$
- $B$ 的质心： $\overset{ˉ}{b} = (\frac{0 + 0 + 1}{3}, \frac{0 + 1 + 0}{3}, \frac{0 + 0 + 0}{3}) = (\frac{1}{3}, \frac{1}{3}, 0)$
去中心化：
- $A$ 去中心化后： ${(- \frac{1}{3}, - \frac{1}{3}, 0), (\frac{2}{3}, - \frac{1}{3}, 0), (- \frac{1}{3}, \frac{2}{3}, 0)}$
- $B$ 去中心化后： ${(- \frac{1}{3}, - \frac{1}{3}, 0), (- \frac{1}{3}, \frac{2}{3}, 0), (\frac{2}{3}, - \frac{1}{3}, 0)}$
计算协方差矩阵 $H$ ：
- $H = \sum_{i} (a_{i} - \overset{a}{ˉ}) (b_{i} - \overset{ˉ}{b})^{T}$
- 对应点对： $A_{1} \leftrightarrow B_{1}$ ， $A_{2} \leftrightarrow B_{3}$ ， $A_{3} \leftrightarrow B_{2}$
- 计算每个点对的外积并求和，得到 $H = \frac{2}{9} - \frac{2}{9} 0 - \frac{2}{9} \frac{2}{9} 0 000$
SVD分解：
- 对 $H$ 进行 SVD 分解，得到 $U$ 和 $V$
- $H$ 的特征值和特征向量计算后， $R = V U^{T} = 010 - 1 00 001$ （绕 z 轴旋转 90 度）
计算平移向量：
- $t = \overset{ˉ}{b} - R \overset{a}{ˉ}$
- $R \overset{a}{ˉ} = 010 - 1 00 001 \frac{1}{3} \frac{1}{3} 0 = - \frac{1}{3} \frac{1}{3} 0$
- $t = (\frac{1}{3}, \frac{1}{3}, 0) - (- \frac{1}{3}, \frac{1}{3}, 0) = (\frac{2}{3}, 0, 0)$
结果：
- 最优旋转矩阵 $R = 010 - 1 00 001$
- 最优平移向量 $t = (\frac{2}{3}, 0, 0)$

验证：应用变换后， $A$ 的点将旋转并平移到与 $B$ 接近重合的位置（由于点对匹配，此处为近似解）。

4. 常见考点与易错点

最近点匹配的高效实现（如k-d树）
SVD分解步骤与旋转矩阵构造
ICP收敛性与初始位置敏感性
多组点云全局配准与误差累积
英文术语拼写与公式记忆

5. 英文关键词与表达

3D data registration, point cloud, rigid transformation, rotation, translation, iterative closest point (ICP), singular value decomposition (SVD), covariance matrix, alignment, convergence

如需更详细例题推导或某一算法的代码实现，可随时补充！

CE7453 Numerical Algorithms Notes