Mar 10, 2026

miniyuan

误差分析

---

题 1.1

题中所给数均加星号表示估计值。

数	绝对误差限	相对误差限
$a^*=3580$	$\vert e(a^*) \vert = 0.5$	$\vert e_r(a^*) \vert = \vert e(a^*)/a^* \vert = 1.397 \times 10^{-4}$
$b^*=0.00476$	$\vert e(b^*) \vert = 5 \times 10^{-6}$	$\vert e_r(b^*) \vert = \vert e(b^*)/b^* \vert = 1.050 \times 10^{-3}$
$c^*=2958 \times 10^{-2}$	$\vert e(c^*) \vert = 5 \times 10^{-3}$	$\vert e_r(c^*) \vert = \vert e(c^*)/c^* \vert = 1.690 \times 10^{-4}$
$d^*=0.1430 \times 10^8$	$\vert e(d^*) \vert = 5 \times 10^{3}$	$\vert e_r(d^*) \vert = \vert e(d^*)/d^* \vert = 3.497 \times 10^{-4}$

题 1.2

证明：

令 $y(x) = \sqrt{x}$，则由误差传播：

$$e_r(y^*) = \frac{dy}{dx}(x^*) \cdot \frac{x^*}{y^*} \cdot e_r(x^*) = \frac{1}{2} \sqrt{x^*} \cdot \frac{x^*}{\sqrt{x^*}} \cdot e_r(x^*) = \frac{1}{2} e_r(x^*)$$

题 1.3

记正方形边长为 $x$，则面积为 $S(x) = x^2$

面积的绝对误差限：

$$\vert e(S^*) \vert = \vert \frac{dS}{dx}(x^*) \cdot e(x^*) \vert = 2 x^* \vert e(x^*) \vert$$

故允许的边长测量绝对误差限为：

$$\vert e(x^*) \vert = \frac{\vert e(S^*) \vert}{2 x^*} = 5 \times 10^{-3} [cm]$$

题 1.4

距离的绝对误差限：

$$\vert e(s^*) \vert = \vert \frac{ds}{dt}(t^*) \cdot e(t^*) \vert = g t^* \vert e(t^*) \vert$$

$e(t^*)$ 不变，$t^*$ 增加时，距离的绝对误差限也增加。

距离的相对误差限：

$$\vert e_r(s^*) \vert = \vert \frac{ds}{dt}(t^*) \cdot \frac{e(t^*)}{s^*} \vert = \frac{2 \vert e(t^*) \vert}{t^*}$$

$e(t^*)$ 不变，$t^*$ 增加时，距离的相对误差限减少。

题 1.5

记 $a = \sqrt{783}$，$a^* = 27.983$，则

$$\vert e(a*) \vert \le \frac{1}{2} \times 10^{-3}$$ $$x_{1, 2} = \frac{56 \pm \sqrt{56^2 - 4}}{2} = 28 \pm a$$

故

$$x_1^* = 28 + a^* = 55.983$$

且绝对误差限

$$\vert e(x_1^*) \vert = \vert e(a^*) \vert \le \frac{1}{2} \times 10^{-3}$$

故 $x_1^*$ 至少有5位有效数字

另一方面，由 $x_1 \cdot x_2 = 1$得：

$$x_2^* = \frac{1}{x_1^*} = 0.017863$$

且绝对误差限

$$\vert e(x_2^*) \vert = \vert -\frac{1}{(x_1^*)^2} \cdot e(x_1^*) \vert \le 1.60 \times 10^{-7}$$

故 $x_2^*$ 至少有5位有效数字

题 1.6

问题 1

我们先对连续世界定义两种圆度：等周型圆度与曲率型圆度，并比较它们之间的差异。本质上等周型圆度是利用曲线的整体性质（周长、面积）来度量，而曲率型圆度则是利用曲线的局部性质（曲率）来度量。所以等周型圆度对震荡不敏感，而后者对曲线的震荡更加敏感。

设 $\gamma \subset \mathbb{R}^2$ 为简单闭曲线（$C^2$ 光滑），并取弧长参数 $s \in [0,L]$。定义：

• $L = \int_\gamma ds$：周长
• $A = \frac{1}{2}\oint_\gamma (x\,dy - y\,dx)$：包围面积
• $\kappa(s)$：相对曲率

定义一：等周型圆度（Isoperimetric Circularity）

$$\xi_{\text{iso}}(\gamma) := \frac{L^2}{4A}$$

定理 1.1（等周不等式）：对任意简单闭曲线 $\gamma$，

$$\xi_{\text{iso}}(\gamma) \geq \pi$$

等号成立当且仅当 $\gamma$ 为圆。

证明：由等周不等式 $L^2 \geq 4\pi A$ 可得。∎

定义二：曲率型圆度（Curvatural Circularity）

设曲率半径 $\rho(s) := \kappa(s)^{-1}$（约定 $\kappa \neq 0$ 对简单闭曲线），定义：

$$\xi_{\text{cur}}(\gamma) := \pi \cdot \frac{\displaystyle\int_0^L \rho(s)^2\,ds}{\left(\displaystyle\int_0^L \rho(s) \,ds\right)^2} \cdot L = \pi \cdot \frac{\mathbb{E}[\rho^2]}{\mathbb{E}[\rho]^2}$$

定理 1.2：对任意简单闭曲线 $\gamma$（$\kappa \neq 0$），

$$\xi_{\text{cur}}(\gamma) \geq \pi$$

等号成立当且仅当 $\rho(s) \equiv \text{const}$，即 $\gamma$ 为圆。

证明：由柯西-施瓦茨不等式，

$$\left(\int \rho\right)^2 \leq \left(\int 1^2\right)\left(\int \rho^2\right) = L \int \rho^2$$

故 $\frac{L \cdot \int \rho^2}{(\int \rho)^2} \geq 1$。等号成立当且仅当 $\rho(s) \equiv c$。∎

近圆凸曲线下两种定义的比较

既然给出了两种定义，那么我们还是希望比较一下两者的差异。利用一点变分法的思想，在近圆曲线下比较一下两种圆度与 $\pi$ 的差异。如果不想看过于复杂的数学推导，可以直接跳到比较部分

设 $\gamma_\epsilon$ 为近圆曲线： $h(\theta) = r + \epsilon u(\theta), \quad \theta \in [0, 2\pi]$

其中：

• $r > 0$ 为半径（常数）
• $\epsilon$ 为小量（常数）
• $u(\theta)$ 为光滑扰动函数，满足标准化条件： $\int_0^{2\pi} u(\theta)\,d\theta = 0$

对凸闭曲线，各几何量可用 $h(\theta)$ 表示为：

• 曲率半径：$\rho(\theta) = h(\theta) + h''(\theta)$
• 周长：$L = \int_0^{2\pi} h(\theta)\,d\theta$
• 面积：$A = \frac{1}{2}\int_0^{2\pi} (h^2 - h'^2)\,d\theta$

1. 等周型圆度

近圆曲线周长：

$$L = \int_0^{2\pi} (r + \epsilon u)\,d\theta = 2\pi r + \epsilon \int_0^{2\pi} u\,d\theta = 2\pi r$$

又因为：

$$h^2 = r^2 + 2\epsilon r u + \epsilon^2 u^2$$ $$h'^2 = \epsilon^2 u'^2$$ $$h^2 - h'^2 = r^2 + 2\epsilon r u + \epsilon^2(u^2 - u'^2)$$

故面积：

$$A = \frac{1}{2}\int_0^{2\pi}\left(r^2 + 2\epsilon r u + \epsilon^2(u^2 - u'^2)\right)d\theta$$ $$= \pi r^2 + \epsilon r \underbrace{\int_0^{2\pi} u\,d\theta}_{=0} + \frac{\epsilon^2}{2}\int_0^{2\pi}(u^2 - u'^2)\,d\theta$$

从而等周圆度：

$$\xi_{\text{iso}} = \frac{L^2}{4A} = \frac{4\pi^2 r^2}{4\left(\pi r^2 + \frac{\epsilon^2}{2}\int(u^2-u'^2)d\theta \right)}$$ $$= \pi \cdot \left(1 + \frac{\epsilon^2}{2\pi r^2}\int_0^{2\pi}(u^2-u'^2)\,d\theta\right)^{-1}$$

对 $\epsilon \ll 1$ 展开也即：

$$\xi_{\text{iso}} = \pi\left(1 - \frac{\epsilon^2}{2\pi r^2}\int_0^{2\pi}(u^2-u'^2)\,d\theta + O(\epsilon^4)\right)$$

1. 曲率型圆度

曲率半径：

$$\rho = h + h'' = r + \epsilon(u + u'') =: r + \epsilon v$$

其中 $v := u + u''$，且 $\int_0^{2\pi} v\,d\theta = \int_0^{2\pi}(u + u'')\,d\theta = 0$（因 $u$ 是周期的）。

计算 $\int_0^L \rho\,ds$（分母）：

$$\int_0^L \rho\,ds = \int_0^{2\pi} \rho^2\,d\theta$$ $$= \int_0^{2\pi}(r + \epsilon v)^2\,d\theta$$ $$= 2\pi r^2 + 2\epsilon r \underbrace{\int_0^{2\pi} v\,d\theta}_{=0} + \epsilon^2 \int_0^{2\pi} v^2\,d\theta$$ $$= 2\pi r^2 + \epsilon^2 \int_0^{2\pi} v^2\,d\theta$$

计算 $\int_0^L \rho^2\,ds$（分子）：

$$\int_0^L \rho^2\,ds = \int_0^{2\pi} \rho^3\,d\theta$$ $$= \int_0^{2\pi}(r + \epsilon v)^3\,d\theta$$ $$= \int_0^{2\pi}(r^3 + 3\epsilon r^2 v + 3\epsilon^2 r v^2 + \epsilon^3 v^3)\,d\theta$$ $$= 2\pi r^3 + 3\epsilon r^2 \underbrace{\int_0^{2\pi} v\,d\theta}_{=0} + 3\epsilon^2 r \int_0^{2\pi} v^2\,d\theta + O(\epsilon^3)$$ $$= 2\pi r^3 + 3\epsilon^2 r \int_0^{2\pi} v^2\,d\theta + O(\epsilon^3)$$

计算 $v^2$ 的积分：

$$\int_0^{2\pi} v^2\,d\theta = \int_0^{2\pi}(u^2 - 2u'^2 + u''^2)\,d\theta$$

其中用到分部积分 $\int uu'' = -\int u'^2$

故曲率圆度：

$$\xi_{\text{cur}} = \pi \cdot \frac{L \cdot \int_0^L \rho^2\,ds}{\left(\int_0^L \rho\,ds\right)^2} = \pi \cdot \frac{2\pi r \cdot \left(2\pi r^3 + 3\epsilon^2 r \int v^2\right)}{\left(2\pi r^2 + \epsilon^2 \int v^2\right)^2}$$ $$= \pi \cdot \frac{4\pi^2 r^4\left(1 + \frac{3\epsilon^2}{2\pi r^2}\int v^2\right)} {4\pi^2 r^4\left(1 + \frac{\epsilon^2}{\pi r^2}\int v^2 + O(\epsilon^4)\right)}$$ $$= \pi \cdot \left(1 + \frac{3\epsilon^2}{2\pi r^2}\int v^2\right)\left(1 - \frac{\epsilon^2}{\pi r^2}\int v^2 + O(\epsilon^4)\right)$$ $$= \pi \cdot \left(1 + \frac{3\epsilon^2}{2\pi r^2}\int v^2 - \frac{2\epsilon^2}{2\pi r^2}\int v^2 + O(\epsilon^4)\right)$$ $$= \pi \cdot \left(1 + \frac{\epsilon^2}{2\pi r^2}\int_0^{2\pi} v^2\,d\theta + O(\epsilon^4)\right)$$

也即：

$$\xi_{\text{cur}} = \pi\left(1 + \frac{\epsilon^2}{2\pi r^2}\int_0^{2\pi}(u^2 - 2u'^2 + u''^2)\,d\theta + O(\epsilon^4)\right)$$

1. 比较

$$\xi_{\text{iso}} = \pi\left(1 - \frac{\epsilon^2}{2\pi r^2}\int_0^{2\pi}(u^2-u'^2)\,d\theta + O(\epsilon^4)\right)$$ $$\xi_{\text{cur}} = \pi\left(1 + \frac{\epsilon^2}{2\pi r^2}\int_0^{2\pi}(u^2 - 2u'^2 + u''^2)\,d\theta + O(\epsilon^4)\right)$$

可以看出等周型圆度中仅包含 $u$ 和 $u'$，而曲率型圆度还包含 $u''$，故前者对曲线的震荡并不敏感，而后者对曲线的震荡更加敏感。事实上，这也是符合我们的直观认识的，对于一个”锯齿形”的近圆曲线，它的周长和面积受到影响并不大，但是曲率就会有很大影响。本质上是因为周长和面积是曲线的整体性质，而曲率则是曲线的局部性质。

---

问题 2

想要将连续情形转化为离散情形，我们需要先定义离散世界中的两种圆度。

设 $S \subset \mathbb{Z}^2$ 为单连通像素区域，$\partial S$ 为其边界像素集。 将 $\partial S$ 中的像素按逆时针排序得到有序边界点序列 $\{p_0, p_1, ..., p_{n-1}\}$，其中 $p_i = (x_i, y_i) \in \mathbb{Z}^2$。

基本几何量定义

面积 $A$：

$$A = \frac{1}{2} | \sum_{i=0}^{n - 1} (x_i y_{i+1} - x_{i+1} y_i) |$$

其中 $(x_{n}, y_{n}) = (x_0, y_0)$

def _polygon_area(verts: List[Tuple[int, int]]) -> float:
    """多边形面积"""
    s = 0
    n = len(verts)
    for i in range(n):
        x1, y1 = verts[i]
        x2, y2 = verts[(i+1) % n]
        s += x1*y2 - x2*y1 # verts 为逆时针点列，故为正
    return abs(s) / 2.0

周长 $P$：

边界点间的欧氏距离之和：

$$P = \sum_{i=0}^{n-1} \|p_{i+1} - p_i\|, \quad p_n = p_0$$

def _exact_perimeter(verts: List[Tuple[int, int]]) -> float:
  """精确周长"""
  per = 0.0
  n = len(verts)
  for i in range(n):
      x1, y1 = verts[i]
      x2, y2 = verts[(i+1) % n]
      per += math.hypot(x2-x1, y2-y1)
  return per

质心 $c = (c_x, c_y)$：

$$\begin{aligned} c_x &= \frac{1}{6A}\sum_{i=0}^{n-1} (x_i + x_{i+1})(x_i y_{i+1} - x_{i+1} y_i) \\ c_y &= \frac{1}{6A}\sum_{i=0}^{n-1} (y_i + y_{i+1})(x_i y_{i+1} - x_{i+1} y_i) \end{aligned}$$

其中 $(x_{n}, y_{n}) = (x_0, y_0)$

def _polygon_centroid(verts: List[Tuple[int, int]]) -> Tuple[float, float]:
  n = len(verts)
  if n == 0:
      return (0, 0)

  area = 0
  cx = cy = 0

  for i in range(n):
      x0, y0 = verts[i]
      x1, y1 = verts[(i+1) % n]
      cross = x0*y1 - x1*y0
      area += cross
      cx += (x0 + x1) * cross
      cy += (y0 + y1) * cross

  if area == 0:
      # 退化多边形，返回顶点平均值
      xs = [v[0] for v in verts]
      ys = [v[1] for v in verts]
      return (sum(xs)/n, sum(ys)/n)

  area *= 0.5
  cx /= (6 * area)
  cy /= (6 * area)
  return (cx, cy)

离散等周型圆度

为了消除局部锯齿效应对周长的影响，我们使用步长 $k$ 采样：

$$\xi_{\text{iso}}^{\text{digital}}(S, k) := \frac{P_k^2}{4A_k}$$

其中：

• $P_k = \frac{1}{k}\sum_{j=0}^{k-1} \sum_{i=0}^{m_j-1} \|p_{i+1}^{(j)} - p_i^{(j)}\|$（多个起始位置平均）
• $A_k = \frac{1}{k}\sum_{j=0}^{k-1} \text{Area}(\{p_i^{(j)}\})$（对应采样点的多边形面积）
• $\{p_i^{(j)}\}$ 是从 $p_j$ 开始每隔 $k$ 个点采样的子序列

def isoperimetric_circularity(self) -> float:
  """
  使用最优步长计算等周型圆度
  """
  if self.optimal_stride is None:
      self.find_optimal_stride()

  if self.optimal_stride == 1:
      return (self.perimeter ** 2) / (4.0 * self.area)

  return self.optimal_xi

最优步长：

$$k_{\text{opt}}(S) = \arg\min_{k \ge 1} |\xi_{\text{iso}}^{\text{digital}}(S, k) - \pi|$$

在实际计算中测试得出最优步长用以计算圆度。

def find_optimal_stride(self, max_stride: int = None) -> Tuple[int, float]:
  """
  寻找最优步长：使圆度最接近π的步长
  返回: (最优步长, 对应的圆度)
  """
  verts = self.ordered
  n = len(verts)

  best_stride = 1
  best_xi = (self.perimeter ** 2) / (4.0 * self.area) # 初始圆度
  best_error = abs(best_xi - math.pi)

  # 测试不同步长
  for stride in range(2, max_stride + 1):
      if n < stride + 2:
          continue

      total_xi = 0.0
      count = 0

      # 多个起始位置取平均
      for start in range(stride):
          sampled = []
          for i in range(start, n, stride):
              sampled.append(verts[i])

          if len(sampled) < 3:
              continue

          # 计算采样多边形的周长和面积
          per = 0.0
          area = 0.0
          m = len(sampled)

          for i in range(m):
              x1, y1 = sampled[i]
              x2, y2 = sampled[(i+1) % m]
              per += math.hypot(x2 - x1, y2 - y1)
              area += x1 * y2 - x2 * y1

          area = abs(area) / 2.0
          if area > 0:
              xi = (per * per) / (4.0 * area)
              total_xi += xi
              count += 1

      if count > 0:
          avg_xi = total_xi / count
          error = abs(avg_xi - math.pi)

          if error < best_error:
              best_error = error
              best_xi = avg_xi
              best_stride = stride

  self.optimal_stride = best_stride
  self.optimal_xi = best_xi
  return best_stride, best_xi

离散曲率型圆度

$$\xi_{\text{cur}}^{\text{digital}}(S) := \pi \cdot \frac{P \cdot \displaystyle\sum_{i=0}^{n-1} d(p_i)^2 \cdot \Delta s_i} {\left(\displaystyle\sum_{i=0}^{n-1} d(p_i) \cdot \Delta s_i\right)^2}$$

其中：

• $d(p_i) = \|p_i - c\|$：边界点到质心的距离
• $\Delta s_i = \|p_{i+1} - p_i\|$：相邻边界点间的弧长微元

def curvature_circularity(self) -> float:
  """曲率型圆度"""
  verts = self.ordered
  n = len(verts)
  if n < 3:
      return float('inf')

  cx, cy = self.centroid
  sum_rho_ds = 0.0
  sum_rho2_ds = 0.0

  for i in range(n):
      x, y = verts[i]
      rho = math.hypot(x - cx, y - cy)
      x2, y2 = verts[(i+1) % n]
      ds = math.hypot(x2 - x, y2 - y)

      sum_rho_ds += rho * ds
      sum_rho2_ds += (rho * rho) * ds

  if sum_rho_ds == 0:
      return float('inf')

  return math.pi * self.perimeter * sum_rho2_ds / (sum_rho_ds ** 2)

最小值

实际最小值应当是与像素网格边长有关的量。由于想不到精确的数学公式，故采用上述离散圆度定义进行数值计算，拟合最小值与误差函数曲线。完整代码见 GitHub 仓库，下面是运行结果：

也即在 $[-R, R] \times [-R, R]$ 的离散像素网格中，拟合得到的等周型圆度最小值：

$$\xi_{\text{iso}}^{min} = \pi + \frac{0.671942}{R^{0.9968}}$$

拟合得到的曲率型圆度最小值：

$$\xi_{\text{cur}}^{min} = \pi + \frac{0.208273}{R^{1.9950}}$$

最后的结果竟然是等周型圆度的误差符合 $\frac{1}{R}$ 变化规律，曲率型圆度的误差符合 $\frac{1}{R^2}$ 变化规律，这其中是否有什么奥秘？还是留待以后探索吧。