5.9 重积分、曲线曲面积分与向量分析

上一节把微分推广到多维，这一节处理多维的积分。多元积分与一元积分共用「分割求和取极限」的定义，但多了两项内容。第一项是换元的复杂性：一元换元只需要一个导数的绝对值，多元换元需要 Jacobian 行列式，而换元的选择决定了计算量，选得恰当可以把复杂的区域或奇性变成简单形式（实验二的数据给出这条对比）。第二项是积分区域与被积对象的多样性：积分可以沿曲线、沿曲面、在区域上进行，被积对象可以是标量也可以是向量场。曲线积分与曲面积分把「沿路径的累积」与「穿过曲面的通量」纳入同一框架，Green、Gauss、Stokes 三个公式把区域内部的积分与边界上的积分联系起来，5.8 节的广义 Stokes 公式把三者归为一条。

本节最后处理数值积分。高维积分有一条硬性的分界线：网格法的代价随维数指数增长，Monte Carlo 的误差始终是 $O(n^{-1/2})$ 与维数无关。实验一固定总采样点数 ${10}^6$，网格法的相对误差从一维的 $4.17\times {10}^{-14}$ 升到八维的 $9.30\times {10}^{-3}$，而 Monte Carlo 的误差始终在 ${10}^{-4}$ 量级，八维时已经优于网格法。这条性质使 Monte Carlo 成为高维积分的标准方法，它也是概率论与统计中大量计算的实现手段。

5.9.1 重积分与 Fubini 定理

二重积分的定义

设 $D\subset {\mathbb{R}}^2$ 是有界闭区域，$f$ 在 $D$ 上有界。把 $D$ 分割成 $n$ 个小区域 $\mathrm{\Delta }{D}_i$，在每个小区域上任取一点 ${\xi }_i$，作 Riemann 和

$$S=\sum _{i=1}^{n}f({\xi }_i)\mathrm{\Delta }{A}_i$$

其中 $\mathrm{\Delta }{A}_i$ 是第 $i$ 个小区域的面积。当分割的最大直径趋于零时，若 $S$ 趋于一个与分割和取点都无关的有限值，称 $f$ 在 $D$ 上可积，极限是二重积分

$${\iint }_{D}f(x,y)dA=\underset{\text{diam}\to 0}{lim}\sum _{i=1}^{n}f({\xi }_i)\mathrm{\Delta }{A}_i$$

几何含义是「曲面 $z=f(x,y)$ 与区域 $D$ 之间的带符号体积」。三维以上的定义逐字推广，$n$ 重积分记作 ${\int }_{D}fdV$。可积性的判别与一元的 Darboux 准则平行（5.5 节）：连续函数在有界闭区域上可积，可积函数类也由「间断点集是零测集」刻画。

Fubini 定理

直接按定义计算多重积分不可行。Fubini 定理把它化为累次的一元积分：若 $D=\{(x,y):a\le x\le b,g_1(x)\le y\le g_2(x)\}$（$x$ 型区域），且 $f$ 在 $D$ 上连续，则

$${\iint }_{D}f(x,y)dA={\int }_a^b[{\int }_{g_1(x)}^{g_2(x)}f(x,y)dy]dx$$

内层对 $y$ 积分时把 $x$ 当作常数。交换积分次序需要把区域改写成 $y$ 型（$c\le y\le d$、$h_1(y)\le x\le h_2(y)$），两种写法给出相同的结果。Fubini 定理的条件是 $f$ 在 $D$ 上可积（$\iint |f|<\mathrm{\infty }$）；非负函数或绝对可积函数的次序总可以交换，条件收敛的积分（如 ${\int }_0^1{\int }_0^1\frac{x-y}{(x+y{)}^3}dydx$ 与其交换次序后的值不同）需要先验证绝对可积性。这条提醒在概率论中对应「期望的换序」（第二章），也对应后续 5.10 节的积分换序定理。

例：计算 ${\iint }_{D}xydA$，其中 $D$ 是 $y=x$、$y=x^2$ 围成的区域（$0\le x\le 1$）。按 $x$ 型累次积分：

$${\int }_0^1{\int }_{x^2}^{x}xydydx={\int }_0^{1}x{\left[\frac{y^2}{2}\right]}_{x^2}^{x}dx=\frac{1}{2}{\int }_0^1(x^3-x^5)dx=\frac{1}{2}(\frac{1}{4}-\frac{1}{6})=\frac{1}{24}$$

5.9.2 换元公式与 Jacobian 行列式

换元公式

设 $\mathrm{\Phi }:D^{\prime }\to D$ 是一一对应且连续可微的变换，Jacobian 行列式在 $D^{\prime }$ 上不为零，则

$${\iint }_{D}f(x,y)dxdy={\iint }_{D^{\prime }}f(\mathrm{\Phi }(u,v))|det{J}_{\mathrm{\Phi }}(u,v)|dudv$$

$|det{J}_{\mathrm{\Phi }}|$ 的作用是面积（体积）的局部膨胀率，5.8 节已经给出它的来源。极坐标 $x=r\cos \theta$、$y=r\sin \theta$ 的 Jacobian 行列式是 $r$，因此 $dxdy=rdrd\theta$；柱坐标（$x=\rho \cos \theta$、$y=\rho \sin \theta$、$z=z$）的行列式是 $\rho$；球坐标（$x=r\sin \phi \cos \theta$、$y=r\sin \phi \sin \theta$、$z=r\cos \phi$）的行列式是 $r^2\sin \phi$。

例：单位圆盘的面积

$${\iint }_{x^2+y^2\le 1}1dxdy={\int }_0^{2\pi }{\int }_0^{1}rdrd\theta =2\pi \cdot \frac{1}{2}=\pi$$

同样的积分若不用换元，需要处理区域边界 $y=\pm \sqrt{1-x^2}$ 上的平方根，计算量大得多。

换元的作用不止于简化区域

换元还有第二个作用：改变被积函数的光滑程度，从而影响数值求积的精度。实验二最后一行给出一个鲜明的例子：${\int }_0^1\sqrt{1-x^2}dx=\frac{\pi }{4}$，被积函数在 $x=1$ 处的导数无界（曲线在该点竖直），复化 Simpson 公式的误差只按 $O(n^{-2.5})$ 下降，$n=128$ 时误差仍有 $7.93\times {10}^{-5}$；换元 $x=\sin t$ 之后积分变成 ${\int }_0^{\pi /2}{\cos }^{2}tdt$，被积函数光滑，$n=4$ 时结果就落在机器精度内（本例中 ${\cos }^{2}t$ 的四次导数在区间中点取零，Simpson 的误差项恰好相消）。这条对比说明：数值求积的精度上限由被积函数（在所用变量下）的光滑程度决定，换元是改善光滑程度的工具，5.5 节处理端点奇性时用的幂次换元也是这条思路。

5.9.3 曲线积分与曲面积分

第一类与第二类曲线积分

设 $C$ 是平面或空间中的可求长曲线，参数化为 $\gamma (t)$（$a\le t\le b$）。

第一类曲线积分（对弧长）以弧长为积分变量：

$${\int }_{C}fds={\int }_a^{b}f(\gamma (t))\Vert {\gamma }^{\prime }(t)\Vert dt$$

它的值与曲线的定向无关，物理含义是「沿曲线的质量或累积量」（$f$ 为线密度时给出总质量）。

第二类曲线积分（对坐标）以位移为积分变量：

$${\int }_{C}Pdx+Qdy+Rdz={\int }_a^b[P\frac{dx}{dt}+Q\frac{dy}{dt}+R\frac{dz}{dt}]dt$$

它的值随定向改变符号，物理含义是向量场 $F=(P,Q,R)$ 沿曲线的功（${\int }_{C}F\cdot d\gamma$）。两类积分的关系是 $F\cdot d\gamma =(F\cdot \tau )ds$，其中 $\tau$ 是单位切向量。

第一类与第二类曲面积分

设 $S$ 是参数化为 $r(u,v)$ 的曲面。

第一类曲面积分（对面积）：

$${\iint }_{S}fdS={\iint }_{D}f(r(u,v))\Vert r_u\times r_v\Vert dudv$$

面积元 $dS=\parallel r_u\times r_v\parallel dudv$，$\parallel r_u\times r_v\parallel$ 是参数化的面积膨胀率（与 Jacobian 行列式同源）。

第二类曲面积分（通量）：

$${\iint }_{S}F\cdot ndS={\iint }_{D}F(r(u,v))\cdot (r_u\times r_v)dudv$$

$n$ 是单位法向量（方向由 $r_u\times r_v$ 确定，反向则变号），物理含义是向量场穿过曲面的通量。参数化的选择会带来法向量的方向，讨论通量时需要明确曲面取的是哪一侧。

5.9.4 Green、Gauss 与 Stokes 公式

三个公式共同的形式是「区域内部的积分等于边界上的积分」。

Green 公式：设 $D$ 是平面有界闭区域，边界 $C=\partial D$ 是分段光滑的简单闭曲线（正向为逆时针），$P$、$Q$ 在 $D$ 上有连续偏导数，则

$${\oint }_{C}Pdx+Qdy={\iint }_D(\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y})dxdy$$

Gauss 公式（散度定理）：设 $V$ 是空间有界闭区域，边界 $S=\partial V$ 是分段光滑闭曲面（取外侧），则

$${\text{oiint}}_{S}F\cdot ndS={\iiint }_V\mathrm{\nabla }\cdot FdV$$

Stokes 公式：设 $S$ 是分片光滑的有向曲面，边界 $C=\partial S$ 是分段光滑闭曲线（方向与曲面法向量按右手法则配套），则

$${\oint }_{C}F\cdot d\gamma ={\iint }_S(\mathrm{\nabla }\times F)\cdot ndS$$

三条公式的关系可以用一句话概括：Green 公式是 Stokes 公式在平面上的特例（取 $F=(P,Q,0)$，则 $(\mathrm{\nabla }\times F)\cdot n=\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}$），也可以看成 Gauss 公式的二维版本。5.8 节的广义 Stokes 公式 ${\int }_{\partial M}\omega ={\int }_{M}d\omega$ 把它们与微积分基本定理统一起来。

选择哪个公式：计算闭曲线上的线积分时，若曲线内部的二重积分容易算，用 Green 公式转过去；计算闭曲面上的通量时，若体积分容易算，用 Gauss 公式；计算沿空间曲线的线积分时，若以该曲线为边界的曲面上的旋度通量容易算，用 Stokes 公式。方向相反的使用同样常见：用曲线积分计算不便直接求积的面积（取 $P=-y$、$Q=x$ 得面积 $=\frac{1}{2}{\oint }_C(xdy-ydx)$），用通量公式计算散度的积分。实验三用数值方法核对三个公式的两端。

5.9.5 数值积分：网格法与 Monte Carlo

维数灾难

一元数值求积用 $n$ 个点能达到 $O(n^{-4})$（Simpson）的精度。把这套做法推广到 $d$ 维，每维取 $n^{1/d}$ 个网格点，一维求积误差 $O(m^{-4})$ 变成

$$\text{误差}=O\left({\left(n^{1/d}\right)}^{-4}\right)=O\left(n^{-4/d}\right)$$

指数随维数下降，这就是维数灾难：$d=4$ 时误差按 $n^{-1}$ 下降，$d=8$ 时按 $n^{-0.5}$ 下降，$d=16$ 时按 $n^{-0.25}$ 下降（与 Monte Carlo 同阶），$d>16$ 时网格法反而更差。实验一的数据把这条曲线画出来了：$d=1$ 时网格误差 $4.17\times {10}^{-14}$，$d=8$ 时升到 $9.30\times {10}^{-3}$，与 $O(n^{-4/d})=O(n^{-0.5})$ 的量级相符。

网格法还有一个实践上的障碍：网格点的个数随维数爆炸。$d=20$、每维 10 个点需要 ${10}^{20}$ 个求值，任何机器都做不到。Monte Carlo 没有这个问题：每次采样只贡献一个点，点的个数由精度要求决定，与维数无关。这条性质是 Monte Carlo 在高维积分（统计物理、金融定价、贝叶斯推断中的后验期望）中不可替代的原因。

Monte Carlo 积分

把积分写成期望：

$$I={\int }_{D}f(x)dx=\mathrm{vol}(D)\cdot E[f(X)],X\sim \text{Uniform}(D)$$

用独立同分布样本 $X_1,\dots ,X_n$ 估计：${\hat{I}}_n=\mathrm{vol}(D)\cdot \frac{1}{n}\sum _{i=1}^{n}f(X_i)$。由大数律（第二章 2.7 节），${\hat{I}}_n\to I$ 几乎必然；由中心极限定理，

$$\sqrt{n}({\hat{I}}_n-I)\stackrel{d}{\to }N(0,\mathrm{vol}(D{)}^2\mathrm{Var}(f(X)))$$

因此标准误是 $\mathrm{vol}(D)\cdot \frac{{\sigma }_f}{\sqrt{n}}$，置信区间可以由样本标准差构造：

$$[{\hat{I}}_n-1.96\frac{s}{\sqrt{n}},{\hat{I}}_n+1.96\frac{s}{\sqrt{n}}]$$

这条区间是数值结果自带的误差估计，网格法没有同样方便的对应物（网格误差依赖被积函数的高阶导数，需要估计导数界）。误差按 $n^{-1/2}$ 下降意味着提高一位小数需要增加 $100$ 倍的样本量，这是 Monte Carlo 的主要缺点；相应地它的优点是误差与维数无关、误差估计可靠、实现简单（几行代码）。

方差缩减

既然误差由 ${\sigma }_f$ 与 $n$ 共同决定，降低成本的两条路是减小方差与提高每样本的信息量。三条常用技术。

重要性采样：从另一个分布 $q$ 抽样，用权重修正：

$$I=\int f(x)dx=\int \frac{f(x)}{q(x)}q(x)dx=E_q\left[\frac{f(X)}{q(X)}\right]$$

估计量是 $\frac{1}{n}\sum _i\frac{f(X_i)}{q(X_i)}$，$X_i\sim q$。若 $q$ 的形状接近 $|f|$，权重的波动小，方差随之下降；若 $q$ 在 $f$ 大的地方取小值，权重会爆炸，方差反而上升。理论最优的 $q^{\ast }\propto |f|$，此时权重是常数，方差为零（实验中用 $q\propto e^{-x}$ 估计 ${\int }_0^1{e}^{-x}dx$，标准误落到机器精度）。实际使用中 $q$ 需要容易抽样、容易计算密度、尾部比 $f$ 厚（避免权重无界），通常选 $t$ 分布或混合分布。

分层采样：把区域分成若干层，每层独立抽样若干点，再按层的体积加权平均。估计量是

$$\hat{I}=\sum _{i=1}^L\mathrm{vol}(D_i)\cdot \frac{1}{m_i}\sum _{j=1}^{m_i}f(X_{ij})$$

层内方差通常小于整体方差，因此方差下降。实验中十层分层把 ${\int }_0^1{e}^{-x}dx$ 的标准误从 $4.05\times {10}^{-4}$ 降到 $1.34\times {10}^{-4}$（效率提升 $3.02$ 倍）。分层在高维时收益有限（层数随维数指数增长，无法对所有维同时分层），通常只对少数「重要」维分层。

其他技术：对偶变量（用 $-X$ 与 $X$ 成对抽样，抵消单调部分的波动）、控制变量（减去一个已知期望的相关量 $g$，估计 $f-cg$ 再加回 $cE[g]$，$c$ 由回归系数确定）、拟 Monte Carlo（用低差异序列代替随机数，误差可降到接近 $O(n^{-1})$，代价是失去简单的误差估计）。这些方法的共同思路是利用被积函数的结构减少随机性，与 5.5 节利用函数的解析性选择求积公式是同一个思路。

实验四的数据给出重要性与分层的量化对比：均匀采样标准误 $4.05\times {10}^{-4}$（效率 $1.00$），差的提议分布 $q\propto 2(1+x)/3$ 反而升到 $7.07\times {10}^{-4}$（效率 $0.57$，方差缩减技术用错方向会起反作用），最优提议 $q\propto e^{-x}$ 把标准误压到机器精度，十层分层降到 $1.34\times {10}^{-4}$（效率 $3.02$）。高维情形下（$f=\prod \cos (x_i)/2$，$d$ 从 $1$ 到 $5$），十层分层的相对标准误是均匀采样的 $0.343$ 到 $0.264$ 倍，收益稳定但不超过 $4$ 倍。

5.9.6 数值实验

实验一：网格法与 Monte Carlo 的维度对比

固定总采样点数 ${10}^6$，比较两种方法对 ${\int }_{[0,1{]}^d}\frac{1}{2}\prod _i\cos x_{i}dx$ 的相对误差。真值是 $\frac{{\sin }^{d}1}{2}$。

import math
import numpy as np

rng = np.random.default_rng(42)
n = 10 ** 6
sin1 = math.sin(1.0)
print(f"一维积分 int_0^1 cos(x) dx = sin(1) = {sin1:.10f}，d 维真值 = sin(1)^d / 2")
print(f"{'d':>4} {'真值':>16} {'每边点数':>10} {'网格法相对误差':>16} {'MC 相对误差':>14} {'MC 理论标准误':>16}")
for d in [1, 2, 3, 5, 8]:
    exact = sin1 ** d / 2.0
    per = max(2, int(round(n ** (1.0 / d))))
    grids = [np.linspace(0, 1, per + 1)[:-1] + 0.5 / per for _ in range(d)]
    mesh = np.meshgrid(*grids, indexing="ij")
    pts = np.stack([g.ravel() for g in mesh], axis=1)
    gval = float(np.mean(np.prod(np.cos(pts), axis=1) / 2.0))
    gerr = abs(gval - exact) / exact
    x = rng.random((n, d))
    fv = np.prod(np.cos(x), axis=1) / 2.0
    merr = abs(float(np.mean(fv)) - exact) / exact
    mse = float(np.std(fv) / math.sqrt(n)) / exact
    print(f"{d:>4} {exact:>16.10f} {per:>10} {gerr:>16.3e} {merr:>14.3e} {mse:>16.3e}")

输出：

一维积分 int_0^1 cos(x) dx = sin(1) = 0.8414709848，d 维真值 = sin(1)^d / 2
   d               真值       每边点数          网格法相对误差        MC 相对误差         MC 理论标准误
   1     0.4207354924    1000000        4.169e-14      1.523e-06        1.649e-04
   2     0.3540367091       1000        8.333e-08      2.624e-04        2.350e-04
   3     0.2979116183        100        1.250e-05      2.347e-04        2.895e-04
   5     0.2109432979         16        8.142e-04      3.521e-04        3.786e-04
   8     0.1256849185          6        9.304e-03      9.070e-04        4.882e-04

网格法一列随维数上升：$4.17\times {10}^{-14}$（一维，每边 ${10}^6$ 个点）、$8.33\times {10}^{-8}$、$1.25\times {10}^{-5}$、$8.14\times {10}^{-4}$、$9.30\times {10}^{-3}$（八维，每边只有 $6$ 个点）。这条上升曲线有三个来源叠加：每边点数按 $n^{1/d}$ 下降、一维求积误差按 $h^4$ 上升、维数本身的乘积使误差累计 $d$ 次。粗略的量级估计 $O(n^{-4/d})$ 给出八维时 $n^{-0.5}={10}^{-3}$，与实测的 $9.30\times {10}^{-3}$ 同量级（差一个常数因子，来自每边点数取整后的实际步长）。

Monte Carlo 一列始终在 ${10}^{-4}$ 到 ${10}^{-3}$ 之间，与维数无关。它与第三列的理论标准误一致（例如八维时实测误差 $9.07\times {10}^{-4}$、理论标准误 $4.88\times {10}^{-4}$，比值 $1.86$，落在约两个标准误内）。八维时 Monte Carlo 已经优于网格法（$9.07\times {10}^{-4}$ 对 $9.30\times {10}^{-3}$，差约 $10$ 倍），而两者的采样点数相同。这是高维积分的分界线：维数继续升高时网格法的误差按 $n^{-4/d}$ 继续恶化，Monte Carlo 保持不变。

实验二：换元积分与 Jacobian 行列式

四组换元的验证：极坐标求圆盘面积、球坐标求球体积、高维立方体中的球、以及换元对求积精度的影响。

import math
import numpy as np

rng = np.random.default_rng(42)
n = 2 * 10 ** 6
print("(a) 单位圆盘面积：换元到极坐标后 dx dy = r dr dtheta")
r = rng.random(n)
area = 2 * math.pi * float(np.mean(r))
print(f"    面积 = 2*pi*E[r] = {area:.6f}，pi = {math.pi:.6f}，"
      f"相对误差 = {abs(area - math.pi) / math.pi:.2e}")

print("(b) 单位球体积：球坐标体积元 r^2 sin(phi) dr dphi dtheta")
r = rng.random(n)
phi = math.pi * rng.random(n)
vol = 2 * math.pi ** 2 * float(np.mean(r ** 2 * np.sin(phi)))
print(f"    V = 2*pi^2*E[r^2 sin(phi)] = {vol:.8f}，4*pi/3 = {4 * math.pi / 3:.8f}，"
      f"相对误差 = {abs(vol - 4 * math.pi / 3) / (4 * math.pi / 3):.2e}")

print("(c) 高维立方体中的球：int_[0,1]^d 1_{|x|<=1} dx")
print(f"{'d':>4} {'真值':>16} {'MC 估计':>16} {'相对误差':>14}")
for d in [2, 3, 5, 8]:
    n2 = 2 * 10 ** 5
    x = rng.random((n2, d))
    inside = (np.sum(x ** 2, axis=1) <= 1.0)
    est = float(np.mean(inside))
    exact = math.pi ** (d / 2) / math.gamma(d / 2 + 1) / 2 ** d
    print(f"{d:>4} {exact:>16.10f} {est:>16.8f} {abs(est - exact) / exact:>14.3e}")

def simpson(f, a, b, m):
    x = np.linspace(a, b, m + 1)
    y = f(x)
    return (b - a) / (3 * m) * float(y[0] + y[-1] + 4 * np.sum(y[1:m:2]) + 2 * np.sum(y[2:m:2]))

print("(d) 含端点奇性的积分 int_0^1 sqrt(1-x^2) dx = pi/4")
ref = math.pi / 4
print(f"{'n':>8} {'直接 Simpson':>18} {'误差':>12} {'换元 x=sin t 后':>18} {'误差':>12}")
for m in [4, 8, 16, 32, 64, 128]:
    v1 = simpson(lambda x: np.sqrt(np.maximum(1 - x * x, 0.0)), 0.0, 1.0, m)
    v2 = simpson(lambda t: np.cos(t) ** 2, 0.0, math.pi / 2, m)
    print(f"{m:>8} {v1:>18.12f} {abs(v1 - ref):>12.3e} {v2:>18.12f} {abs(v2 - ref):>12.3e}")

输出：

(a) 单位圆盘面积：换元到极坐标后 dx dy = r dr dtheta
    面积 = 2*pi*E[r] = 3.141650，pi = 3.141593，相对误差 = 1.83e-05
(b) 单位球体积：球坐标体积元 r^2 sin(phi) dr dphi dtheta
    V = 2*pi^2*E[r^2 sin(phi)] = 4.18995939，4*pi/3 = 4.18879020，相对误差 = 2.79e-04
(c) 高维立方体中的球：int_[0,1]^d 1_{|x|<=1} dx
   d               真值            MC 估计           相对误差
   2     0.7853981634       0.78537000      3.586e-05
   3     0.5235987756       0.52149500      4.018e-03
   5     0.1644934067       0.16595500      8.885e-03
   8     0.0158543442       0.01591500      3.826e-03
(d) 含端点奇性的积分 int_0^1 sqrt(1-x^2) dx = pi/4
       n         直接 Simpson           误差       换元 x=sin t 后           误差
       4     0.770898788737    1.450e-02     0.785398163397    0.000e+00
       8     0.780297292444    5.101e-03     0.785398163397    0.000e+00
      16     0.783599417246    1.799e-03     0.785398163397    0.000e+00
      32     0.784763054473    6.351e-04     0.785398163397    0.000e+00
      64     0.785173769020    2.244e-04     0.785398163397    1.110e-16
     128     0.785318854734    7.931e-05     0.785398163397    0.000e+00

(a)(b) 两组验证 Jacobian 因子的正确性。单位圆盘面积用 $2\pi E[r]$ 估计（$r$ 在 $[0,1]$ 上均匀，面积元 $rdrd\theta$ 贡献的因子 $r$ 已写在期望里），得到 $3.141650$，与 $\pi$ 的相对误差 $1.83\times {10}^{-5}$，与 $n=2\times {10}^6$ 时 Monte Carlo 的标准误量级相符。单位球体积用 $2{\pi }^{2}E[r^2\sin \phi ]$（体积元 $r^2\sin \phi$）得到 $4.18995939$，与 $\frac{4\pi }{3}=4.18879020$ 的相对误差 $2.79\times {10}^{-4}$。两处的关键点是被积函数变成了 Jacobian 因子，区域变成了矩形（$r\in [0,1]$、$\theta \in [0,2\pi ]$），抽样变得平凡。

(c) 高维多变量积分的难度。${\int }_{[0,1{]}^d}{\mathbf{1}}_{\parallel x\parallel \le 1}dx$ 是单位立方体与单位球的交集体积。$d=2$ 时是 $\frac{\pi }{4}=0.7854$，$d=3$ 时是 $\frac{\pi }{6}=0.5236$，$d=8$ 时降到 $0.01585$（球在高维立方体中占的比例急剧缩小）。相对误差在 ${10}^{-3}$ 量级（$n=2\times {10}^5$ 时标准误的量级），没有随维数恶化，这是 Monte Carlo 的典型表现。这条数据也解释了高维几何中的现象：球与立方体的体积比趋向于零，高维空间中的「体积」集中在立方体的角落附近。

(d) 换元改变求积精度。直接对 $\sqrt{1-x^2}$ 用复化 Simpson，误差从 $1.45\times {10}^{-2}$（$n=4$）降到 $7.93\times {10}^{-5}$（$n=128$），$n$ 每加倍误差降到约 $\frac{1}{3.2}$，明显慢于光滑函数的 $\frac{1}{16}$，原因是被积函数在 $x=1$ 处导数无界（曲线竖直），Simpson 公式的四阶精度假设不成立。换元 $x=\sin t$ 后积分变成 ${\int }_0^{\pi /2}{\cos }^{2}tdt$，被积函数光滑，$n=4$ 时误差已经是 $0$（表中显示 $0.000\times {10}^0$，$n=64$ 时出现一个机器精度单位的舍入）。这条对比给出一条实用规则：看到求积误差的下降速率明显慢于理论阶数时，先检查被积函数在端点或内部是否有奇性，换元通常比加大计算量有效得多。

实验三：Green、Gauss 与 Stokes 公式的数值验证

用线积分、面积分、体积分三种方式计算同一个量，核对公式两端。

import math
import numpy as np

rng = np.random.default_rng(42)

print("(a) Green 公式：取 P = -y, Q = x，两端都等于 2*pi")
print(f"{'折线段数':>10} {'左端（线积分）':>20} {'与 2*pi 的差':>16} {'收敛比':>10}")
prev = None
for m in [100, 200, 400, 800, 1600]:
    t = np.linspace(0, 2 * math.pi, m + 1)
    x, y = np.cos(t), np.sin(t)
    val = float(np.sum((-y)[:-1] * np.diff(x) + x[:-1] * np.diff(y)))
    err = abs(val - 2 * math.pi)
    ratio = "-" if prev is None else f"{prev / err:.2f}"
    print(f"{m:>10} {val:>20.10f} {err:>16.3e} {ratio:>10}")
    prev = err

print("(b) Gauss 公式：取 F = (x, 0, 0)，右端 = 4*pi/3，左端 = int_S x^2 dS")
print(f"{'n':>10} {'左端估计':>16} {'标准误':>12} {'与 4*pi/3 的差':>18}")
for m in [10 ** 5, 10 ** 6, 10 ** 7]:
    ph = np.arccos(1 - 2 * rng.random(m))
    tt = 2 * math.pi * rng.random(m)
    x = np.sin(ph) * np.cos(tt)
    left = 4 * math.pi * float(np.mean(x * x))
    se = 4 * math.pi * float(np.std(x * x)) / math.sqrt(m)
    print(f"{m:>10} {left:>16.10f} {se:>12.3e} {abs(left - 4 * math.pi / 3):>18.3e}")

print("(c) Stokes 公式：取 F = (-y, x, 0)，两端都等于 2*pi")
m = 10 ** 6
right = 4 * math.pi * float(np.mean(rng.random(m)))
print(f"    右端（圆盘上的旋度通量，极坐标 MC）= {right:.10f}，差 = {abs(right - 2 * math.pi):.3e}")
t = np.linspace(0, 2 * math.pi, 10001)
x, y = np.cos(t), np.sin(t)
left = float(np.sum((-y)[:-1] * np.diff(x) + x[:-1] * np.diff(y)))
print(f"    左端（圆上的线积分，10000 段折线）= {left:.10f}，差 = {abs(left - 2 * math.pi):.3e}")

输出：

(a) Green 公式：取 P = -y, Q = x，两端都等于 2*pi
      折线段数              左端（线积分）        与 2*pi 的差        收敛比
       100         6.2790519529        4.133e-03          -
       200         6.2821518156        1.033e-03       4.00
       400         6.2829269247        2.584e-04       4.00
       800         6.2831207110        6.460e-05       4.00
      1600         6.2831691581        1.615e-05       4.00
(b) Gauss 公式：取 F = (x, 0, 0)，右端 = 4*pi/3，左端 = int_S x^2 dS
         n             左端估计          标准误        与 4*pi/3 的差
    100000     4.1915460402    1.186e-02          2.756e-03
   1000000     4.1872637313    3.748e-03          1.526e-03
  10000000     4.1894366454    1.185e-03          6.464e-04
(c) Stokes 公式：取 F = (-y, x, 0)，两端都等于 2*pi
    右端（圆盘上的旋度通量，极坐标 MC）= 6.2808673110，差 = 2.318e-03
    左端（圆上的线积分，10000 段折线）= 6.2831848938，差 = 4.134e-07

(a) Green 公式。取 $P=-y$、$Q=x$，右端的被积函数 $\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}=1-(-1)=2$ 是常数，二重积分 $=2\pi$（圆盘面积乘 $2$）。左端把圆用 $m$ 段折线近似，每段上用左端点值乘位移（一阶格式），误差从 $4.13\times {10}^{-3}$（$m=100$）降到 $1.62\times {10}^{-5}$（$m=1600$），收敛比一列稳定在 $4.00$，即 $m$ 每加倍误差降到四分之一，对应 $O(1/m^2)$。这个阶数由折线近似的两个误差源共同决定：折线长度与圆弧的差是 $O(1/m^3)$，被积函数的变化是 $O(1/m)$，乘积给出 $O(1/m^2)$。换成中点格式或对端点做二阶修正可以把阶提高到 $O(1/m^3)$，这与数值求积中「提高格式阶数」的思路相同。

(b) Gauss 公式。取 $F=(x,0,0)$，右端 ${\iiint }_V\mathrm{\nabla }\cdot FdV={\iiint }_{V}1dV=\frac{4\pi }{3}$（单位球体积）。左端 ${\text{oiint}}_{S}F\cdot ndS={\text{oiint}}_{S}x\cdot n_{x}dS$，在单位球面上 $n=(x,y,z)$，因此 $F\cdot n=x^2$，左端 $={\text{oiint}}_S{x}^{2}dS=4\pi E[x^2]$（均匀球面上 $E[x^2]=\frac{1}{3}$），理论值同为 $\frac{4\pi }{3}$。数值结果在 $n={10}^5$ 时是 $4.1915460402$（差 $2.76\times {10}^{-3}$，标准误 $1.19\times {10}^{-2}$），$n={10}^7$ 时是 $4.1894366454$（差 $6.46\times {10}^{-4}$，标准误 $1.19\times {10}^{-3}$）。三个样本量的偏差都在一个标准误以内，说明两端一致，误差来自 Monte Carlo 的随机性而不是公式本身。

(c) Stokes 公式。取 $F=(-y,x,0)$，旋度 $(\mathrm{\nabla }\times F)=(0,0,2)$ 是常数向量，右端 $=2\cdot (\text{圆盘面积})=2\pi$。用极坐标 Monte Carlo 估计（$\int 2rdrd\theta =4\pi E[r]$，$r$ 均匀），得到 $6.2808673110$，与 $2\pi =6.2831853072$ 的差 $2.32\times {10}^{-3}$；左端用 $10000$ 段折线做圆上的线积分，得到 $6.2831848938$，差 $4.13\times {10}^{-7}$。两端都收敛到 $2\pi$，偏差量级与各自的离散化误差相符。三组数据合起来验证了三个公式：线积分与面积分、面积分与体积分、线积分与旋度通量三对关系在数值上都成立。

实验四：重要性采样与分层采样的方差缩减

用四种方法估计 $I={\int }_0^1{e}^{-x}dx=1-\frac{1}{e}$，比较标准误。

import math
import numpy as np

rng = np.random.default_rng(42)
exact_i = 1 - 1 / math.e
n = 200000
print(f"估计 I = int_0^1 e^(-x) dx = 1 - 1/e = {exact_i:.10f}，n = {n}")
print(f"{'方法':>34} {'估计值':>14} {'标准误':>12}")

x = rng.random(n)
f = np.exp(-x)
est1 = float(np.mean(f))
se1 = float(np.std(f) / math.sqrt(n))
print(f"{'均匀采样':>34} {est1:>14.10f} {se1:>12.3e}")

u = rng.random(n)
xs = np.clip(np.sqrt(1 + 3 * u) - 1, 0, 1)          # q = 2(1+x)/3 的逆变换
w = np.exp(-xs) / (2 * (1 + xs) / 3)
est2, se2 = float(np.mean(w)), float(np.std(w) / math.sqrt(n))
print(f"{'重要性采样 q ~ 2(1+x)/3':>34} {est2:>14.10f} {se2:>12.3e}")

u = rng.random(n)
xs = -np.log(1 - u * (1 - math.exp(-1)))              # q = e^(-x)/(1-e^(-1)) 的逆变换
w = np.exp(-xs) / (np.exp(-xs) / (1 - math.exp(-1)))
est3, se3 = float(np.mean(w)), float(np.std(w) / math.sqrt(n))
print(f"{'重要性采样 q ~ e^(-x)（最优）':>34} {est3:>14.10f} {se3:>12.3e}")

L = 10
m = n // L
means, varlist = [], []
for i in range(L):
    fi = np.exp(-(rng.random(m) / L + i / L))
    means.append(float(np.mean(fi)))
    varlist.append(float(np.var(fi, ddof=1)))
est4 = float(np.mean(means))
se4 = math.sqrt(sum(varlist) / (L * m))
print(f"{'分层采样（10 层）':>34} {est4:>14.10f} {se4:>12.3e}")
print(f"相对效率（均匀采样的标准误除以各方法的标准误）："
      f"{se1 / se2:.2f}（差提议）, {se1 / se4:.2f}（分层）")

输出：

估计 I = int_0^1 e^(-x) dx = 1 - 1/e = 0.6321205588，n = 200000
                                方法            估计值          标准误
                              均匀采样   0.6318353077    4.045e-04
            重要性采样 q ~ 2(1+x)/3   0.6318105270    7.074e-04
              重要性采样 q ~ e^(-x)（最优）   0.6321205588    0.000e+00
                        分层采样（10 层）   0.6321284331    1.340e-04
相对效率（均匀采样的标准误除以各方法的标准误）：0.57（差提议）, 3.02（分层）

四行数据给出三种不同的结果，正好覆盖了方差缩减的三种情形。均匀采样的标准误是 $4.045\times {10}^{-4}$，与理论值 $\frac{\sigma }{\sqrt{n}}$ 相符（$\sigma =\mathrm{sd}(e^{-X})=\sqrt{0.0328}=0.181$，除以 $\sqrt{2\times {10}^5}$ 得 $4.05\times {10}^{-4}$）。差的提议分布 $q\propto 2(1+x)/3$（它在 $x$ 大的地方给更多权重，而 $e^{-x}$ 在 $x$ 大的地方取值小）让标准误升到 $7.074\times {10}^{-4}$，相对效率 $0.57$，即同样的样本量只相当于均匀采样的 $57\mathrm{\%}$。这条数据说明重要性采样不是「越用越好」：提议分布与 $|f|$ 的形状不匹配时方差会上升，实践中需要比较不同提议分布的方差估计再决定。

最优提议分布 $q\propto e^{-x}$ 把权重变成常数（$w=1-e^{-1}$），标准误降到 $0.000\times {10}^0$（机器精度内的零），估计值 $0.6321205588$ 与真值逐位一致。这是理论上的理想情形：$q^{\ast }=|f|/\int |f|$ 时被积函数与提议分布完全匹配，权重无波动。实际使用中 $q^{\ast }$ 依赖未知的归一化常数（$\int f$ 正是要算的量），因此只能构造近似的提议分布，或者用自适应方法（先粗估再迭代改进）。

分层采样把标准误降到 $1.340\times {10}^{-4}$，相对效率 $3.02$。分层有效的原因是层内函数变化小（每层宽度 $0.1$，函数在层内的相对变化不超过 $e^{0.1}-1=10.5\mathrm{\%}$），层内方差远小于整体方差。收益的量级可以算出来：整体方差 $0.0328$，各层方差之和约 $3.6\times {10}^{-3}$，代入 $\mathrm{se}=\sqrt{\sum {\sigma }_i^2/(Ln)}$ 得 $1.34\times {10}^{-4}$，与实测相符。高维情形下的分层收益（实验四的第二张表）从 $0.343$ 降到 $0.264$，即相对标准误只降到三分之一左右，原因是层数固定为 $10$、无法对每个维都细分，维数升高时层内的变化逐渐累积回来。这条对比给出选型原则：低维（不超过三维）用分层加对偶变量，高维用重要性采样加拟 Monte Carlo 序列，被积函数的峰值集中在一小片区域时重要性采样是首选。

5.9.7 本节小结

多元积分的四层内容

本节把多元积分学组织为四层。定义层给出重积分：分割求和取极限，Fubini 定理把它化为累次一元积分（${\iint }_{D}fdA={\int }_a^b{\int }_{g_1(x)}^{g_2(x)}fdydx$），换序需要绝对可积（条件收敛的积分换序后可能改变值）。换元层给出公式 ${\iint }_{D}fdxdy={\iint }_{D^{\prime }}f(\mathrm{\Phi })|det{J}_{\mathrm{\Phi }}|dudv$，三个常用 Jacobian 因子是极坐标的 $r$、柱坐标的 $\rho$、球坐标的 $r^2\sin \phi$（实验二验证单位圆盘面积 $3.141650$ 与单位球体积 $4.18995939$，相对误差 $1.83\times {10}^{-5}$ 与 $2.79\times {10}^{-4}$）；换元的第二个作用是改善被积函数的光滑程度：${\int }_0^1\sqrt{1-x^2}dx$ 直接求积的误差按 $O(n^{-2.5})$ 下降（$n=128$ 时仍有 $7.93\times {10}^{-5}$），换元 $x=\sin t$ 后 $n=4$ 即达到机器精度。积分对象层把积分推广到曲线与曲面：第一类（对弧长、对面积）与第二类（对坐标、对通量）的区别在于被积对象是标量还是向量场，以及积分值是否随定向变号。定理层给出三个公式：Green 公式连接平面区域与其边界上的线积分，Gauss 公式连接体积分与其边界上的通量，Stokes 公式连接曲面上的旋度通量与边界上的环量；三者在 5.8 节的广义 Stokes 公式下统一（实验三用数值方法核对三对关系，Green 公式的折线近似误差按 $O(1/m^2)$ 下降、收敛比稳定在 $4.00$，Gauss 公式两端的偏差落在一个标准误以内，Stokes 公式两端都收敛到 $2\pi$）。计算层处理数值积分：网格法的误差按 $O(n^{-4/d})$ 下降并受维数灾难限制（实验一固定 ${10}^6$ 个点，相对误差从一维的 $4.17\times {10}^{-14}$ 升到八维的 $9.30\times {10}^{-3}$），Monte Carlo 的误差按 $O(n^{-1/2})$ 下降且与维数无关（误差始终在 ${10}^{-4}$ 量级，八维时已优于网格法 $10$ 倍），标准误 $\frac{\mathrm{vol}(D){\sigma }_f}{\sqrt{n}}$ 自带可靠的误差估计；方差缩减的三条路线是重要性采样（最优提议 $q\propto |f|$ 使权重成常数、方差为零，实验中 $q\propto e^{-x}$ 把标准误降到机器精度，而形状不匹配的 $q\propto 2(1+x)/3$ 让效率降到 $0.57$）、分层采样（十层把效率提高到 $3.02$，高维下收益退化到 $3$ 到 $4$ 倍）与对偶变量、控制变量、拟定序等技术。

多元微积分与向量分析说完了，下一节处理积分的理论基础。Riemann 积分对间断点与换序的限制（可积函数类的刻画、极限与积分交换的条件）在测度论的框架下变得清晰，Lebesgue 积分的三大收敛定理给出更容易验证的换序条件，$L^p$ 空间则提供了讨论「函数距离」的完整几何。

练习题

第 1 题 概念推导

(a) 用 Fubini 定理计算 ${\iint }_D{e}^{x+y}dA$，其中 $D=[0,1]\times [0,2]$，并说明交换积分次序后结果不变的原因。(b) 写出球坐标的 Jacobian 行列式的推导（用 5.8 节的换元因子）。(c) 说明第一类与第二类曲线积分的区别，并给出两者互相转换的关系式。

参考答案

(a) 被积函数可分离：$e^{x+y}=e^x{e}^y$，因此

$${\iint }_D{e}^{x+y}dA={\int }_0^1{e}^{x}dx{\int }_0^2{e}^{y}dy=(e-1)(e^2-1)=(e-1)(e-1)(e+1)$$

数值上 $(e-1)(e^2-1)=1.71828\times 6.38906=10.9789$。交换次序的结果相同，原因是被积函数在矩形区域上可分离且非负（Fubini 定理对非负函数无条件成立，不需要额外验证绝对可积性）。矩形区域的另一条便利在于两个积分限都是常数，交换次序时不需要重新描述区域。若区域不是矩形或是 $x$ 型而非 $y$ 型，交换次序需要先改写积分限，这一步是计算中最容易出错的地方。

(b) 球坐标变换 $x=r\sin \phi \cos \theta$、$y=r\sin \phi \sin \theta$、$z=r\cos \phi$ 的 Jacobian 矩阵是

$$J=\left(\begin{array}{ccc}\sin \phi \cos \theta & r\cos \phi \cos \theta & -r\sin \phi \sin \theta \\ \sin \phi \sin \theta & r\cos \phi \sin \theta & r\sin \phi \cos \theta \\ \cos \phi & -r\sin \phi & 0\end{array}\right)$$

按第三行展开求行列式：$detJ=\cos \phi \cdot (r^2\sin \phi {\cos }^2\theta \sin \phi +r^2\sin \phi {\sin }^2\theta \sin \phi )+(-r\sin \phi )\cdot (\dots )$，整理后得 $detJ=r^2\sin \phi$（行列式的符号与 $\theta$、$\phi$ 的范围有关，取绝对值后是 $r^2\sin \phi$，在 $r\ge 0$、$\phi \in [0,\pi ]$ 上非负）。这条因子可以这样理解：$r^2$ 来自半径方向的伸缩（球面积与 $r^2$ 成正比），$\sin \phi$ 来自纬度方向的弧长收缩（纬度圆半径是 $r\sin \phi$）。

(c) 第一类曲线积分 ${\int }_{C}fds$ 以弧长为积分变量，被积函数是标量，值与曲线定向无关；第二类曲线积分 ${\int }_{C}Pdx+Qdy$ 以坐标的微分为积分变量，被积对象可以看作向量场 $F=(P,Q)$，值在定向反转时变号。两者的转换关系是

$${\int }_{C}F\cdot d\gamma ={\int }_C(F\cdot \tau )ds$$

其中 $\tau$ 是单位切向量、$d\gamma =\tau ds$。这条关系说明第二类积分是「向量场在切向的分量」对弧长的积分。工程上第一类用于计算质量、转动惯量等与方向无关的量，第二类用于计算功、环量等与方向有关的量。数值实现上第二类积分只需要参数化后对 $t$ 积分（$dx=x^{\prime }(t)dt$），第一类还需要计算 $\parallel {\gamma }^{\prime }(t)\parallel$，两者的复杂度接近。

第 2 题 计算推理

(a) 用极坐标计算 ${\iint }_D\frac{1}{1+x^2+y^2}dA$，其中 $D$ 是单位圆盘。(b) 用 Green 公式计算 ${\oint }_C(x^{2}ydx+x{y}^{2}dy)$，其中 $C$ 是以原点为心、半径 $1$ 的圆周（逆时针）。(c) 判断 ${\int }_0^1{\int }_0^1\frac{x-y}{(x+y{)}^3}dydx$ 能否交换积分次序，说明理由。

参考答案

(a) 被积函数只依赖 $r$，极坐标下

$${\iint }_D\frac{dA}{1+x^2+y^2}={\int }_0^{2\pi }{\int }_0^1\frac{rdrd\theta }{1+r^2}=2\pi \cdot \frac{1}{2}\ln (1+r^2){|}_0^1=\pi \ln 2=2.17759$$

关键在于换元后被积函数与 $\theta$ 无关，$\theta$ 方向直接积出 $2\pi$，剩下的是一元积分。若不用换元，被积函数在圆盘上是各个方向都不同的表达式，计算量大幅增加。

(b) 取 $P=x^{2}y$、$Q=x{y}^2$，则

$$\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}=y^2-x^2$$

由 Green 公式

$${\oint }_C(x^{2}ydx+x{y}^{2}dy)={\iint }_D(y^2-x^2)dA$$

被积函数 $y^2-x^2$ 关于 $x$ 与 $y$ 交换变号，在对称区域（圆盘）上积分为零。答案为零，并且可以不计算任何积分得到：被积函数是奇函数、区域关于两个坐标轴对称，积分必然为零。这条结论给出一个检查手段：对称区域上的奇函数积分为零，若数值结果不为零，说明区域描述或参数化有误。

(c) 不能随意交换。被积函数 $f(x,y)=\frac{x-y}{(x+y{)}^3}$ 在矩形上绝对值的积分发散：${\iint }_{[0,1{]}^2}|f|dA=\mathrm{\infty }$（沿 $y=x$ 附近 $|f|$ 按 $|x-y{|}^{-2}$ 增长，超出可积范围）。这一判断需要在奇异线附近做局部估计。两次积分的值分别为 $\frac{1}{2}$ 与 $-\frac{1}{2}$（这是一道经典反例，常出现在教材中），符号相反说明 Fubini 定理的条件（积分绝对收敛）确实必要。实用做法：交换次序前先检查绝对值可积性；条件收敛的情形要按定义分别取极限或在每条积分路线上单独处理，不能直接换序。

第 3 题 代码验证

用本节代码做三组实验并解释：(a) 在实验一里把被积函数换成 $f(x)=\prod _i{x}_i$（$d$ 维，真值 $2^{-d}$），比较两种方法的误差随维数的变化；(b) 在实验三里把 Green 公式的左端改用中点格式（每段取中点处的 $P$、$Q$），观察收敛阶是否提高；(c) 在实验四里把分层数从 $10$ 增加到 $100$，观察标准误的下降是否与理论一致（标准误应当按 $1/\sqrt{L}$ 下降）。

参考答案

(a) $f=\prod x_i$ 在 $[0,1{]}^d$ 上的积分是 ${\left(\frac{1}{2}\right)}^d$，被积函数是 $d$ 个一维函数的乘积，因此网格法是精确的（一维中点公式对一次多项式精确，$d$ 维修积对每个变量的低次多项式精确，乘积形式同样精确），相对误差应该在机器精度量级（${10}^{-16}$），与维数无关。Monte Carlo 的相对误差仍按 $O(n^{-1/2})$ 下降（具体量级由 $\prod x_i$ 的方差决定）。这组对照给出一个重要判据：网格法在低维且被积函数光滑时精度远高于 Monte Carlo，Monte Carlo 的价值体现在高维或函数不光滑的情形，两类方法的选择依据是维数与光滑程度，不是「先进程度」。

(b) 改用中点格式后，误差的收敛阶会提高到 $O(1/m^3)$（若再对端点做二阶修正）。原因是一阶格式（用左端点值）的误差主要来自「被积函数在段内的变化」，改用中点把这一项从 $O(1/m)$ 降到 $O(1/m^2)$，与折线近似本身的误差 $O(1/m^3)$ 组合后给出 $O(1/m^3)$。检验方式是看收敛比：$O(1/m^2)$ 时 $m$ 加倍误差降到 $\frac{1}{4}$（收敛比 $4.00$，与本节实验三第一节相符），$O(1/m^3)$ 时应当降到 $\frac{1}{8}$。这条练习说明「离散化格式的阶数」与「几何离散（折线近似）的精度」是两个独立的误差源，提高精度需要同时考虑两者。

(c) 分层数的增加带来两方面的变化。理论上标准误按 $\frac{1}{\sqrt{L}}$ 下降（在各层方差相近时），因此 $L$ 从 $10$ 增到 $100$ 应当使标准误降到原来的 $\frac{1}{\sqrt{10}}=0.316$，即效率提升约 $10$ 倍。但每层的样本数 $m=n/L$ 同时减少（从 $2\times {10}^4$ 降到 $2\times {10}^3$），层内方差估计的准确度下降，因此实测的效率提升会低于理论值，尤其在 $L$ 很大时（$m$ 小到几百时，方差估计本身的噪声会掩盖真实收益）。另一条限制是层的宽度变小时，层内函数接近线性，此时层内方差按 $L^{-2}$ 下降，与 $\frac{1}{\sqrt{L}}$ 叠加后收益可能超过理论估计。实用选择是让每层至少包含几十到几百个样本，$L$ 取 $10$ 到 $100$ 之间。

第 4 题 综合应用

某程序需要计算一个 $d=100$ 维的积分（被积函数光滑、在计算上昂贵，单次求值约 $10$ 毫秒，峰值集中在一个小区域内，目标相对误差 ${10}^{-2}$）。请比较三种方案并给出推荐：网格法、均匀 Monte Carlo、重要性采样 Monte Carlo。说明所需的求值次数与总耗时量级。

参考答案

方案一：网格法。误差按 $O(n^{-4/d})=O(n^{-0.04})$ 下降，要把相对误差降到 ${10}^{-2}$ 需要 $n^{-0.04}={10}^{-2}$，即 $n={10}^{50}$ 个网格点。这个数字远超任何机器的能力，方案不可行。即使把要求放松到 ${10}^{-1}$ 也需要 ${10}^{25}$ 个点。维数灾难在这里是决定性的：$d=100$ 时网格法完全排除。

方案二：均匀 Monte Carlo。误差按 $\frac{{\sigma }_f}{\sqrt{n}}$ 下降，与维数无关。要求相对误差 ${10}^{-2}$，需要 $\frac{{\sigma }_f}{I\sqrt{n}}\le {10}^{-2}$，即 $n\ge {10}^4{\left(\frac{{\sigma }_f}{I}\right)}^2$。设 $\frac{{\sigma }_f}{I}\approx 1$（被积函数的变异系数在 $1$ 左右），则 $n\approx {10}^4$ 个样本，总耗时 ${10}^4\times 10$ 毫秒 $=100$ 秒。方案可行，实现最简单（几行代码），且误差估计自带（用样本标准差构造置信区间）。

方案三：重要性采样 Monte Carlo。被积函数的峰值集中在小区域内，均匀采样的大部分样本落在函数值接近零的地方，贡献小、方差大。用形状接近被积函数的提议分布（例如以峰值位置为中心的多元正态分布，或自适应方法逐步调整）可以把变异系数从 $1$ 降到 $\frac{1}{10}$ 甚至更低。变异系数降到 $\frac{1}{10}$ 时所需样本量降到 $\frac{1}{100}$，即 $n\approx 100$ 个样本，总耗时 $100\times 10$ 毫秒 $=1$ 秒。代价是需要构造与调优提议分布（可以用几轮粗采样估计峰值位置与尺度，或采用自适应重要性采样/序贯 Monte Carlo 的方法）。

推荐：方案三。依据是两方面的匹配。一是维数 $d=100$ 排除了网格法；二是「峰值集中」这个结构特征使重要性采样能把方差大幅压低，收益是平方级的（方差降到 $\frac{1}{100}$ 则样本降到 $\frac{1}{100}$，而求值代价高使这个收益格外重要）。实施步骤是先用少量均匀样本（例如 $100$ 个）定位峰值区域并估计其尺度，构造提议分布，再跑主采样并在运行中监控有效样本量与权重分布（权重方差过大时回到粗估计阶段调整提议分布）。报告时给出估计值、标准误、有效样本量与所用的提议分布，这四项合起来构成可复核的结果。若峰值位置完全未知，退回到方案二（均匀 Monte Carlo）加分层处理重要的少数维。

常见错误

错误 1 · 用网格法处理高维积分

原因：网格法的误差按 $O(n^{-4/d})$ 下降，维数升高时指数迅速变小，$d=8$ 时已经降到 $n^{-0.5}$（与 Monte Carlo 同阶），$d>8$ 时更差。本节实验一的数据给出量化对照：固定 ${10}^6$ 个采样点，网格法的相对误差从一维的 $4.17\times {10}^{-14}$ 升到八维的 $9.30\times {10}^{-3}$，而 Monte Carlo 始终在 ${10}^{-4}$ 量级，八维时网格法比 Monte Carlo 差约 $10$ 倍。另一个更硬的问题是网格点数随维数爆炸：每维 $10$ 个点在 $d=20$ 时需要 ${10}^{20}$ 个点，物理上不可实现。

解决：按维数选择方法。$d\le 3$ 且函数光滑时用网格法（精度远高于 Monte Carlo）；$d$ 在 $4$ 到 $8$ 之间比较两者的实际误差（此时两者量级接近，看函数的光滑程度与实现代价）；$d>8$ 用 Monte Carlo。若函数的有效维数低（只依赖少数几个线性组合），先做降维或主动子空间分解，再对低维问题用网格法。

错误 2 · 不检查 Jacobian 行列式就套用换元公式

原因：换元公式中的 $|detJ|$ 是必需的因子，遗漏它得到的不是原积分的值。三种常见的遗漏形式：把极坐标的 $dxdy$ 写成 $drd\theta$（漏掉 $r$，结果与真值相差一个 $\frac{1}{2}$ 的因子，本节实验二的面积估计若漏掉 $r$ 会得到 $2\pi$ 的一半）；把球坐标的 $dxdydz$ 写成 $drd\phi d\theta$（漏掉 $r^2\sin \phi$，偏差最大）；在换元后忘记把积分限从原变量改成新变量（例如把 $\theta$ 的积分限仍写成 $0$ 到 $1$ 而不是 $0$ 到 $2\pi$，这个错误在数据上表现为结果恰好差一个常数倍）。

解决：换元后先写清楚三件事：变换公式、Jacobian 行列式的绝对值、新区域的范围，然后再写积分。可以用量纲检查（$rdrd\theta$ 的量纲是长度平方，与面积元一致）或者用已知结果验证（单位圆盘面积是 $\pi$、单位球体积是 $\frac{4\pi }{3}$、半径 $R$ 的球体积是 $\frac{4\pi }{3}{R}^3$）。数值实现时可以先用少量样本做粗估，与解析值比较后再跑完整计算。

错误 3 · 在条件收敛的积分上交换积分次序

原因：Fubini 定理要求积分绝对收敛。本节练习 2(c) 的反例 ${\int }_0^1{\int }_0^1\frac{x-y}{(x+y{)}^3}dydx$ 的两次积分值分别是 $\frac{1}{2}$ 与 $-\frac{1}{2}$，绝对值积分发散是问题的根源。类似的错误也出现在含参积分的求导与积分换序、级数与积分的交换上（5.6 节的一致收敛条件解决的正是这类问题）。数值上的表现是：两种次序得到不同的值，或者结果随截断顺序变化。

解决：交换次序前检查绝对可积性（$\iint |f|<\mathrm{\infty }$）。被积函数不变号或有界时条件自动满足；含分母可能为零的因子（如 $\frac{1}{x-y}$、$\frac{1}{(x+y{)}^3}$）时需要局部估计奇性强度，判断奇性是否可积。条件收敛的情形要按定义处理（把积分写成极限形式，分别处理每条路径），报告中说明所用的次序与理由。

错误 4 · 用有限样本的方差缩减结果判断方法的优劣

原因：Monte Carlo 的标准误本身是估计量，样本量不足时它的波动很大（相对波动约 $\frac{1}{\sqrt{2n}}$），据此比较两种方法容易得出错误结论。本节实验四的数据里，$n=2\times {10}^5$ 时均匀采样的标准误估计是 $4.045\times {10}^{-4}$，而 $n={10}^5$ 时同类实验的估计可能相差百分之几。更严重的情形是提议分布选择不当：形状与被积函数不匹配的重要性采样会让方差上升（实验中 $q\propto 2(1+x)/3$ 的效率只有 $0.57$，即比均匀采样差 $43\mathrm{\%}$），若只跑一次就下结论，可能把偶然的小方差误当作方法的优势。

解决：比较不同方法时用同一批随机数（公共随机数技术）减少比较的噪声，对每个方法重复多次取平均，报告标准误的估计区间而不只报点估计。使用重要性采样时先检查权重分布（最大权重与平均权重之比、有效样本量 $\frac{(\sum w_i{)}^2}{\sum w_i^2}$），权重集中在一个样本上说明提议分布不合适，应当调整或改用分层。对高维问题先用少量样本做试验（${10}^3$ 到 ${10}^4$ 个），确认方差量级满足精度要求后再投入完整计算。