Apr 11, 2026

miniyuan

解线性方程组的迭代法

---

题 3.1

1. Jacobi 迭代法

   $$\begin{bmatrix} 10 & 0 & 0 \\ 0 & 10 & 0 \\ 0 & 0 & 5 \\ \end{bmatrix} \mathbf{x}^{(k+1)} = \begin{bmatrix} 0 & 2 & 1 \\ 2 & 0 & 1 \\ 1 & 2 & 0 \\ \end{bmatrix} \mathbf{x}^{(k)} + \begin{bmatrix} 3 \\ 15 \\ 10 \\ \end{bmatrix}$$

   也即有迭代公式：

   $$\begin{cases} x_1^{(k+1)} = \dfrac{2x_2^{(k)} + x_3^{(k)} + 3}{10} \\ x_2^{(k+1)} = \dfrac{2x_1^{(k)} + x_3^{(k)} + 15}{10} \\ x_3^{(k+1)} = \dfrac{x_1^{(k)} + 2x_2^{(k)} + 10}{5} \end{cases}$$

   代入 $x^{(0)}=[0,0,0]^\top$ 得 $\mathbf{x}^{(1)} = [3/10,\ 3/2,\ 2]^\top$

   $$\begin{cases} x_1^{(2)} = \frac{2(3/2) + 2 + 3}{10} = \frac{3+2+3}{10} = \frac{4}{5} \\ x_2^{(2)} = \frac{2(3/10) + 2 + 15}{10} = \frac{3/5+2+15}{10} = \frac{44}{25} \\ x_3^{(2)} = \frac{3/10 + 2(3/2) + 10}{5} = \frac{3/10+3+10}{5} = \frac{133}{50} \end{cases}$$

   也即 $\mathbf{x}^{(2)} = [0.8,\ 1.76,\ 2.66]^\top$

2. Gauss-Seidel 迭代法

   $$\begin{bmatrix} 10 & 0 & 0 \\ -2 & 10 & 0 \\ -1 & -2 & 5 \\ \end{bmatrix} \mathbf{x}^{(k+1)} = \begin{bmatrix} 0 & 2 & 1 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \\ \end{bmatrix} \mathbf{x}^{(k)} + \begin{bmatrix} 3 \\ 15 \\ 10 \\ \end{bmatrix}$$

   也即有迭代公式：

   $$\begin{cases} x_1^{(k+1)} = \dfrac{2x_2^{(k)} + x_3^{(k)} + 3}{10} \\ x_2^{(k+1)} = \dfrac{2x_1^{(k+1)} + x_3^{(k)} + 15}{10} \\ x_3^{(k+1)} = \dfrac{x_1^{(k+1)} + 2x_2^{(k+1)} + 10}{5} \end{cases}$$

   代入 $x^{(0)}=[0,0,0]^\top$ 得：

   $$\begin{cases} x_1^{(1)} = \frac{3}{10} \\ x_2^{(1)} = \frac{2(3/10) + 15}{10} = \frac{3/5+15}{10} = \frac{39}{25} \\ x_3^{(1)} = \frac{3/10 + 2(39/25) + 10}{5} = \frac{3/10+78/25+10}{5} = \frac{671}{250} \end{cases}$$ $$\begin{cases} x_1^{(2)} = \frac{2(39/25) + 671/250 + 3}{10} = \frac{2201}{2500} \\ x_2^{(2)} = \frac{2(2201/2500) + 671/250 + 15}{10} = \frac{12153}{6250} \\ x_3^{(2)} = \frac{2201/2500 + 2(12153/6250) + 10}{5} = \frac{184617}{62500} \end{cases}$$

   也即 $\mathbf{x}^{(2)} = [0.8804,\ 1.94448,\ 2.953872]^\top$

题 3.2

迭代格式为：

$$\begin{bmatrix} a_{11} & 0 \\ 0 & a_{22} \end{bmatrix} \mathbf{x}^{(k+1)} = \begin{bmatrix} 0 & -a_{12} \\ -a_{21} & 0 \end{bmatrix} \mathbf{x}^{(k)} + \begin{bmatrix} b_1 \\ b_2 \end{bmatrix}$$

从而迭代矩阵：

$$\mathbf{M} = \begin{bmatrix} 0 & -a_{12}/a_{11} \\-a_{21}/a_{22} & 0 \end{bmatrix}$$

特征根为 $\lambda = \pm \sqrt{\frac{a_{12}a_{21}}{a_{11}a_{22}}}$（注意可能为纯虚数），从而谱半径：

$$\rho(\mathbf{M}) = \sqrt{\left| \frac{a_{12}a_{21}}{a_{11}a_{22}} \right|}$$

迭代收敛的充要条件为 $\sqrt{| \frac{a_{12}a_{21}}{a_{11}a_{22}} |} \lt 1$，也即迭代收敛的充要条件为：

$$\left| \frac{a_{12}a_{21}}{a_{11}a_{22}} \right| \lt 1$$

题 3.3

1. 最速下降法

   编写 python 代码如下：

   import numpy as np
   
   x = np.array([0, 0, 0], dtype=float)
   A = np.array([[4, 3, 0], [3, 4, -1], [0, -1, 4]], dtype=float)
   b = np.array([24, 30, -24], dtype=float)
   r = b - A @ x
   
   print(f"x^(0) = {x}")
   print(f"r^(0) = {r}")
   
   for k in range(3):
       t = A @ r
       alpha = np.dot(r, r) / np.dot(r, t)  # 步长
       x += alpha * r                       # 更新解
       r -= alpha * t                       # 更新残差
       print(f"alpha^({k}) = {alpha:.5f}")
       print(f"x^({k+1}) = {x}")
       print(f"r^({k+1}) = {r}")

   运行结果：

   x^(0) = [0. 0. 0.]
   r^(0) = [ 24.  30. -24.]
   alpha^(0) = 0.14691
   x^(1) = [ 3.5257732   4.40721649 -3.5257732 ]
   r^(1) = [-3.32474227 -1.73195876 -5.48969072]
   alpha^(1) = 0.22980
   x^(2) = [ 2.76173276  4.00920473 -4.78732834]
   r^(2) = [ 0.92545478  0.89065447 -0.84148191]
   alpha^(2) = 0.14852
   x^(3) = [ 2.89917841  4.14148195 -4.9123026 ]
   r^(3) = [-0.02115948 -0.17576563 -0.20930764]

   从而得到

   $$\mathbf{x}^{(3)} \approx [2.8992,\ 4.1415,\ -4.9123]^\top$$

   绝对误差约为

   $$\mathbf{e}^{(3)} \approx [-0.1008,\ 0.1415,\ 0.0877]^\top$$

2. 共轭梯度法

   编写 python 代码如下：

   import numpy as np
   
   x = np.array([0, 0, 0], dtype=float)
   A = np.array([[4, 3, 0], [3, 4, -1], [0, -1, 4]], dtype=float)
   b = np.array([24, 30, -24], dtype=float)
   r = b - A @ x
   d = r.copy()
   r_norm = np.dot(r, r)         # 存储 ||r||^2
   
   print(f"x^(0) = {x}")
   print(f"r^(0) = {r}")
   print(f"d^(0) = {d}")
   
   for k in range(3):
       t = A @ d
       s = np.dot(d, t)
       alpha = np.dot(r, d) / s  # 步长
       x += alpha * d            # 更新解
       r -= alpha * t            # 更新残差
   
       r_norm_new = np.dot(r, r)
       beta = r_norm_new / r_norm
       r_norm = r_norm_new
   
       d = r + beta * d          # 更新搜索方向
       print(f"alpha^({k}) = {alpha:.5f}")
       print(f"beta^({k}) = {beta:.5f}")
       print(f"x^({k+1}) = {x}")
       print(f"r^({k+1}) = {r}")
       print(f"d^({k+1}) = {d}")
   

   运行结果：

   x^(0) = [0. 0. 0.]
   r^(0) = [ 24.  30. -24.]
   d^(0) = [ 24.  30. -24.]
   alpha^(0) = 0.14691
   beta^(0) = 0.02154
   x^(1) = [ 3.5257732   4.40721649 -3.5257732 ]
   r^(1) = [-3.32474227 -1.73195876 -5.48969072]
   d^(1) = [-2.80789669 -1.0859018  -6.0065363 ]
   alpha^(1) = 0.23782
   beta^(1) = 0.00071
   x^(2) = [ 2.85801112  4.14897194 -4.95422216]
   r^(2) = [ 0.1210397  -0.12414328 -0.0341394 ]
   d^(2) = [ 0.11905546 -0.12491065 -0.038384  ]
   alpha^(2) = 1.19263
   beta^(2) = 0.00000
   x^(3) = [ 3.  4. -5.]
   r^(3) = [-1.80411242e-15 -1.44328993e-15 -1.03389519e-15]
   d^(3) = [-1.80411242e-15 -1.44328993e-15 -1.03389519e-15]

   从而得到

   $$\mathbf{x}^{(2)} \approx [2.8580,\ 4.1490,\ -4.9542]^\top$$

   绝对误差约为

   $$\mathbf{e}^{(2)} \approx [-0.1420,\ 0.1490,\ 0.0458]^\top$$

题 3.4

证明：

因为 $\mathbf{A} = \mathbf{D} - \mathbf{L} - \mathbf{U}$，所以：

$$\mathbf{B} = \mathbf{I} - \mathbf{D}^{-1} (\mathbf{D} - \mathbf{L} - \mathbf{U}) = \mathbf{D}^{-1} \mathbf{L} + \mathbf{D}^{-1} \mathbf{U}$$

其中 $\mathbf{D}^{-1} \mathbf{L}$ 为严格下三角矩阵，$\mathbf{D}^{-1} \mathbf{U}$ 为严格上三角矩阵。

于是式 $\text{P3.4-1}$ 可化为：

$$\tilde{\mathbf{x}}^{(k)} = \mathbf{D}^{-1} \mathbf{U} \mathbf{x}^{(k)} + \mathbf{D}^{-1} \mathbf{L} \tilde{\mathbf{x}}^{(k)} + \mathbf{g}$$

也即：

$$\begin{aligned} \tilde{\mathbf{x}}^{(k)} &= (\mathbf{I} - \mathbf{D}^{-1} \mathbf{L})^{-1} (\mathbf{D}^{-1} \mathbf{U} \mathbf{x}^{(k)} + \mathbf{g}) \\ &= (\mathbf{D} - \mathbf{L})^{-1} (\mathbf{U} \mathbf{x}^{(k)} + \mathbf{D} \mathbf{g}) \end{aligned}$$

式 $\text{P3.4-2}$ 可化为：

$$\mathbf{x}^{(k+1)} = \mathbf{D}^{-1} \mathbf{U} \mathbf{x}^{(k+1)} + \mathbf{D}^{-1} \mathbf{L} \tilde{\mathbf{x}}^{(k)} + \mathbf{g}$$

也即：

$$\begin{aligned} \mathbf{x}^{(k+1)} &= (\mathbf{I} - \mathbf{D}^{-1} \mathbf{U})^{-1} (\mathbf{D}^{-1} \mathbf{L} \tilde{\mathbf{x}}^{(k)} + \mathbf{g}) \\ &= (\mathbf{D} - \mathbf{U})^{-1} (\mathbf{L} \tilde{\mathbf{x}}^{(k)} + \mathbf{D} \mathbf{g}) \\ &= (\mathbf{D} - \mathbf{U})^{-1} [\mathbf{L} (\mathbf{D} - \mathbf{L})^{-1} (\mathbf{U} \mathbf{x}^{(k)} + \mathbf{D} \mathbf{g}) + \mathbf{D} \mathbf{g}] \\ &= (\mathbf{D} - \mathbf{U})^{-1} \mathbf{L} (\mathbf{D} - \mathbf{L})^{-1} \mathbf{U} \mathbf{x}^{(k)} + (\mathbf{D} - \mathbf{U})^{-1} [\mathbf{L} (\mathbf{D} - \mathbf{L})^{-1} + \mathbf{I}] \mathbf{D} \mathbf{g} \\ &= (\mathbf{D} - \mathbf{U})^{-1} \mathbf{L} (\mathbf{D} - \mathbf{L})^{-1} \mathbf{U} \mathbf{x}^{(k)} + (\mathbf{D} - \mathbf{U})^{-1} \mathbf{D} (\mathbf{D} - \mathbf{L})^{-1} \mathbf{D} \mathbf{g} \end{aligned}$$

从而

$$\begin{aligned} \mathbf{M} &= (\mathbf{D} - \mathbf{U})^{-1} \mathbf{L} (\mathbf{D} - \mathbf{L})^{-1} \mathbf{U} \\ \mathbf{f} &= (\mathbf{D} - \mathbf{U})^{-1} \mathbf{D} (\mathbf{D} - \mathbf{L})^{-1} \mathbf{D} \end{aligned}$$

证毕。

题 3.5

我们假设己方和己方机器人的任一视线中不存在两个及以上的猎人，否则参考下图可举反例证明无法精确定位。

若某区域恰被 $k$ 个机器人监测，则称其为 $k$ 覆盖的。 我们先证明，在已知猎人数量上界 $n$ 时，我们可以通过 $(2n + 1)$ 覆盖唯一确定所有猎人的位置。 然后再利用这一性质设计确定所有猎人位置的算法，我们首先探测猎人数量的上界，然后再进行精确计算。

引理

已知猎人数量上界为 $n$ 时，我们可以通过 $(2n + 1)$ 覆盖的机器人布置唯一确定所有猎人的位置。 这里要求每个 $(2n + 1)$ 覆盖区域的监测机器人不存在三点共线。

证明：

只需证明对每个 $(2n + 1)$ 覆盖区域中的猎人都可以唯一确定位置， 也即证明 $(2n + 1)$ 覆盖区域中不存在两个不完全相同的猎人集合，使得 $(2n + 1)$ 个监测机器人对他们的视角集合完全相同。

反证法。反设机器人为 $\mathcal{R} = \set{R_1, R_2, \cdots, R_m}$， 两个不同的猎人集合 $\mathcal{P} = \set{P_i, P_2, \cdots, P_n}$ 和 $\mathcal{Q} = \set{Q_i, Q_2, \cdots, Q_n}$， 其中 $m = 2n + 1$。

对每个机器人 $R_k$，存在 $n$ 排列 $\pi_k$ 使得从 $R_k$ 看 $P_i$ 与 $Q_{\pi_k(i)}$ 视角相同，也即 $R_k P_i Q_{\pi_k(i)}$ 三点共线。

定义 $S_{ij} = \{k \in \set{1, \cdots, m} : \pi_k(i) = j\}$，即与 $P_i Q_j$ 共线的机器人集合。

若 $P_i \neq Q_j$，则所有 $S_{ij}$ 中的机器人都在直线 $P_i Q_j$ 上。由机器人中无三点共线得 $|S_{ij}| \leq 2$。

因点集 $\mathcal{P}$ 与点集 $\mathcal{Q}$ 不完全相同，不妨 $P_1 \notin \mathcal{Q}$。

则集合族 $\set{S_{11}, S_{12}, \cdots, S_{1n}}$ 构成了对下标集 $\set{1, \cdots, m}$ 的划分。 也即每个下标 $k \in \set{1, \cdots, m}$ 都恰好属于且仅属于一个集合 $S_{1\pi_k(1)}$。

故有：

$$m = \sum_{j=1}^n |S_{1j}| \leq 2n$$

矛盾！

整体算法

1. 分块： 将整个森林区域分块，划分为若干直径小于 $2d$ 的连通子区域 $D_1, D_2, \cdots, D_s$。
2. 探测猎人上界： 在每个子区域的中心处部署一名机器人并暴露，共计消耗 $s$ 个机器人。区域 $D_i$ 机器人观测到不同方向的数量 $n_i$ 即为该区域猎人数上界。 从而整体上界为 $n = \sum_{i=1}^s n_i$。
3. 精确定位： 对每个子区域周围部署 $m_i \ge 2n_i + 1$ 名机器人，然后再在每个子区域的中心处部署一名机器人并暴露，共计消耗 $s$ 个机器人。由引理知，每个区域可以完成精确定位，从而完成了所有猎人定位。

精确定位算法

对每个非空子区域 $D_i$，已知其猎人数上界 $n_i \leq n$，且周围已部署 $m_i \geq 2n_i+1$ 台无三点共线的观测机器人 $R_1, \ldots, R_{m_i}$。

记第 $k$ 台机器人观测到的无序方向集合为 $\Theta_k = \{\theta_{k,1}, \ldots, \theta_{k,t_k}\}$。

1. 候选点生成：

   对每对 $(\alpha, \beta) \in \Theta_1 \times \Theta_2$，求解：

   $$R_1 + t \begin{pmatrix} \cos\alpha \\ \sin\alpha \end{pmatrix} = R_2 + s \begin{pmatrix} \cos\beta \\ \sin\beta \end{pmatrix}, \quad t, s > 0$$

   这是关于 $(t, s)$ 的 2 阶线性系统：

   $$\begin{pmatrix} \cos\alpha & -\cos\beta \\ \sin\alpha & -\sin\beta \end{pmatrix} \begin{pmatrix} t \\ s \end{pmatrix} = R_2 - R_1$$
   • 若 $\det = \sin(\beta - \alpha) = 0$ 或 $t \leq 0$ 或 $s \leq 0$，舍弃该对。
   • 否则解出 $(t, s)$，计算交点： $$X_{\alpha,\beta} = R_1 + t \begin{pmatrix} \cos\alpha \\ \sin\alpha \end{pmatrix}$$

   对有效交点 $X_{\alpha,\beta}$，验证是否位于所有机器人观测范围内：

   $$\|X_{\alpha,\beta} - R_k\| \leq d, \quad k = 1, 2, \cdots, m_i$$

   剔除不满足距离约束的点。

   进一步验证交点是否位于目标子区域 $D_i$ 内：

   $$X_{\alpha,\beta} \in D_i$$

   最终得到候选点集合：

   $$\mathcal{C} = \left\{ X_{\alpha,\beta} : (\alpha,\beta) \in \Theta_1 \times \Theta_2,\; t,s > 0,\; \|X_{\alpha,\beta}-R_k\| \leq d,\; X_{\alpha,\beta} \in D_i \right\}$$

2. 一致性筛选：

   引入第三台机器人 $R_3$。对每个 $X \in \mathcal{C}$，计算其相对于 $R_3$ 的方位角：

   $$\phi_3(X) = \operatorname{atan2}(X - R_3)$$

   仅保留满足 $\phi_3(X) \in \Theta_3$ 的点，记筛选后集合为 $\mathcal{C}' \subseteq \mathcal{C}$。

   对 $\mathcal{C}'$ 同上利用其余机器人 $R_4, \ldots, R_{m_i}$ 进行筛选，最终得到通过筛选的候选点集合 $\mathcal{P}_i^*$。

3. 唯一性输出：

   由引理知，$\mathcal{P}_i^*$ 恰好唯一地对应真实猎人集合。因此直接输出：

   $$\mathcal{P}_i = \mathcal{P}_i^*$$

区域划分

不妨假设森林边界为矩形 $[0,L]\times[0,W]$，采用正方形网格，边长为 $\sqrt{2} d$，则子区域数量为：

$$s = \left\lceil \frac{L}{\sqrt{2}d} \right\rceil \times \left\lceil \frac{W}{\sqrt{2}d} \right\rceil = O\left(\frac{LW}{d^2}\right)$$

更一般地，最优划分方式为正六边形划分（ 参见 Kershner, 1939 ），此处为便于计算简化考虑。

复杂度

• 分块：$O(s)$
• 探测上界：$O(s \cdot n_{\max})$
• 精确定位：$\displaystyle O\left(\sum_{i=1}^{s} m_i n_i^2\right)$

由引理，精确定位阶段对每个子区域 $D_i$ 部署的观测机器人数量为：

$$m_i = 2n_i + 1$$

其中 $n_i$ 为该区域猎人数上界。因此：

$$\sum_{i=1}^{s} m_i n_i^2 = \sum_{i=1}^{s} (2n_i + 1) n_i^2 = O\left(\sum_{i=1}^{s} n_i^3\right)$$

设森林中猎人总数为 $N$。由 Hölder 不等式：

$$\sum_{i=1}^{s} n_i^3 \leq \left(\max_i n_i\right)^2 \cdot \sum_{i=1}^{s} n_i = n_{\max}^2 N$$

故整体时间复杂度为：

$$O\left(\frac{LW}{d^2} \cdot n_{\max} + n_{\max}^2 N\right)$$

总消耗为：

$$2s = O\left(\frac{LW}{d^2}\right)$$

参考文献

Kershner, R. (1939). The Number of Circles Covering a Set. American Journal of Mathematics, 61(3), 665–671.