高斯积分这个坑，我前前后后踩了不下五次。第一次是在高数下册，盯着 $\int e^{-x^2}dx$ 发呆，老师一句"原函数不是初等函数"就把我们打发了；第二次是在大学物理，麦克斯韦速率分布那个归一化常数 $4\pi(m/2\pi kT)^{3/2}$ 从天而降，没人告诉我它怎么来的；第三次是在概率统计课上，正态分布的密度函数前面那个 $1/\sqrt{2\pi}\sigma$ 又一次凭空出现。直到有一次帮人调试一段数据拟合代码，模型残差始终对不上，最后发现问题出在我把高斯积分的一般形式里的缩放因子写错了——那一刻我才真正把这三门课里的同一个东西串起来。高斯积分，说白了就是形如 $\int_{-\infty}^{+\infty} e^{-ax^2+bx+c}dx$ 这一类积分的统称，它的核心结论只有一个：$\int_{-\infty}^{+\infty}e^{-x^2}dx=\sqrt{\pi}$。就这么一个数，撑起了正态分布、量子谐振子、高斯光束、热统配分函数、卡尔曼滤波、高斯过程回归一大片内容。这篇东西写给三类人：正在被高数折磨、想知道这玩意儿到底有啥用的学生；物理课上被归一化常数搞晕、想弄明白推导链条的同学；以及已经工作、需要在代码里算高斯型积分但总是差个系数的工程实践者。我会把推导、参数、代码、踩坑全摊开讲，能抄作业的地方直接给结论。
1. 一个"积不出来"的积分，逼出了三条完全不同的路
1.1 原函数不是初等函数，这句话到底意味着什么
很多人对"积不出来"这四个字有误解，以为意思是这个积分不存在或者没有值。完全不是。$\int_0^x e^{-t^2}dt$ 对任意有限的 $x$ 都有确定的值，它是一个良好定义的函数，只是这个函数没法用幂函数、指数函数、对数函数、三角函数、反三角函数这五类基本初等函数经过有限次加减乘除和复合表示出来。换句话说，它的原函数存在，但不属于初等函数这个"封闭俱乐部"。
这个事实的证明思路来自微分代数里的刘维尔定理，大意是：如果一个初等函数的原函数仍然是初等函数，那么它必须满足非常苛刻的结构条件。$e^{-x^2}$ 恰好不满足。你在网上看到的所谓"误差函数" $\mathrm{erf}(x)=\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}dt$，本质上就是给这个积不出来的东西起了个名字，然后当成一个新函数用——这不是解决了问题，这是绕过了问题。
我见过不少同学在这里陷入死循环：反复尝试分部积分、换元、有理化，试图把 $e^{-x^2}$ 积出来。浪费时间是小事，更麻烦的是形成一种错误直觉，觉得"只要方法够巧就一定能积出来"。真实情况是，数学里有一大类函数就是这样，你得接受它，然后换工具。接受这件事本身，就是从"计算数学"走向"分析数学"的分水岭。
1.2 三条出路：级数展开、数值逼近、升维打击
既然直着走不通，历史上人们给出了三条完全不同的路线，每一条都对应一种思维方式，值得单独拎出来说。
第一条是级数展开。把 $e^{-x^2}$ 展开成 $\sum_{n=0}^{\infty}(-1)^n x^{2n}/n!$，然后逐项积分，得到 $\sum_{n=0}^{\infty}(-1)^n x^{2n+1}/[n!(2n+1)]$。这条路在 $|x|$ 不大的时候收敛很快，工程上算误差函数基本都用类似的办法。缺点是 $x$ 一大就灾难性抵消，收敛慢得令人绝望，而且它给不出一个漂亮的封闭形式答案。
第二条是数值积分。梯形法、辛普森法、龙贝格积分，或者更聪明的 Gauss 型求积。这条路的优势是通用，什么函数都能算；劣势是它只给你数字，不给你结构。你想知道"为什么归一化常数是 $\sqrt{\pi}$ 而不是 $\sqrt{2\pi}$"，数值积分永远回答不了。
第三条是升维，也就是把一维积分平方成一个二重积分，再换到极坐标下去算。这是最漂亮的一条路，也是高斯积分这个名字真正闪光的时刻。它的精髓在于：一维算不动的东西，放到二维里突然就变成了可以分离变量的形式，而极坐标下的面积元又恰好贡献出一个 $r$，把那个讨厌的平方指数给"吃"掉了。
我第一次看到这个操作时的反应是"这也行？"——因为它本质上是在用几何结构换代数结构。后来做信号处理和机器学习，发现这套思路反复出现：时域卷积难算就换到频域，空域滤波难算就换到小波域。升维不是数学家的炫技，它是一种通用的解题策略。
2. 极坐标那条路：平方、换元、开方，每一步都有讲究
2.1 先把合法性说清楚：二重积分凭什么能拆成两个一维积分的乘积
设 $I=\int_{-\infty}^{+\infty}e^{-x^2}dx$。因为被积函数恒正，$I$ 要么是有限正数，要么发散到无穷，不存在符号问题。现在考虑 $I^2$：
$$I^2=\left(\int_{-\infty}^{+\infty}e^{-x^2}dx\right)\left(\int_{-\infty}^{+\infty}e^{-y^2}dy\right)$$
这里有一个容易糊弄过去但必须交代的步骤：两个积分相乘，能不能写成二重积分 $\iint_{\mathbb{R}^2}e^{-(x^2+y^2)}dxdy$？答案是能，前提是两个积分都绝对收敛。而 $e^{-x^2}$ 在全实轴上的积分确实绝对收敛（被 $e^{-|x|}$ 在 $|x|>1$ 上控制，而 $e^{-|x|}$ 可积），所以富比尼-托内利定理保证了乘积可以写成二重积分。这一步在课堂上经常被一句"显然"带过，但如果你以后要做多重积分的变量替换、或者处理非绝对收敛的情形（比如 $\int \sin x/x$ 那种条件收敛的），这个合法性检查就是救命的东西。
我见过一个非常典型的错误：有人对 $\int_0^\infty \frac{\sin x}{x}dx$ 做平方，然后拆成二重积分换极坐标，结果得出一个错误答案。问题就出在这个积分是条件收敛的，不满足富比尼定理的前提。高斯积分没这个麻烦，但养成"先查绝对收敛"的习惯，能帮你避开后面 80% 的坑。
2.2 雅可比因子 $r$ 不是凭空冒出来的，它是面积的度量
把 $(x,y)$ 换成 $(r,\theta)$，$x=r\cos\theta$，$y=r\sin\theta$，那么 $x^2+y^2=r^2$，被积函数变成 $e^{-r^2}$，干净得让人感动。但面积元不能直接写成 $dr,d\theta$，必须写成 $r,dr,d\theta$。这个 $r$ 是雅可比行列式的绝对值：
$$J=\begin{vmatrix}\frac{\partial x}{\partial r}&\frac{\partial x}{\partial \theta}\\frac{\partial y}{\partial r}&\frac{\partial y}{\partial \theta}\end{vmatrix}=\begin{vmatrix}\cos\theta&-r\sin\theta\\sin\theta&r\cos\theta\end{vmatrix}=r\cos^2\theta+r\sin^2\theta=r$$
从几何上理解更直观：在极坐标里，$r$ 越大，同样的 $d\theta$ 对应的弧长越长（弧长等于 $r,d\theta$），所以面积元自然带着一个 $r$ 的放大因子。你可以想象一个圆环扇区，内半径 $r$、外半径 $r+dr$、张角 $d\theta$，它的面积约等于 $r,dr,d\theta$——当 $r$ 很小的时候这个扇区几乎退化，当 $r$ 很大时它被拉得很开。
这个因子是整条推导的关键。因为 $e^{-r^2}\cdot r$ 恰好凑成了 $-\frac{1}{2}\frac{d}{dr}e^{-r^2}$ 的形式，直接就能积出来。如果没有这个 $r$，你换到极坐标也白搭。我后来做图像处理里的极坐标变换、做雷达信号的距离-方位角重采样，每次都要盯着雅可比因子有没有漏，漏了就是全盘皆错。
2.3 完整链条走一遍，$\sqrt{\pi}$ 是怎么落地的
把上面几步串起来：
$$I^2=\int_{0}^{2\pi}\int_{0}^{+\infty}e^{-r^2}r,dr,d\theta$$
内层对 $r$ 积分：令 $u=r^2$，$du=2r,dr$，于是 $\int_0^{+\infty}e^{-r^2}r,dr=\frac{1}{2}\int_0^{+\infty}e^{-u}du=\frac{1}{2}$。外层对 $\theta$ 积分给出 $2\pi$。所以 $I^2=2\pi\cdot\frac{1}{2}=\pi$，即 $I=\sqrt{\pi}$。
注意最后开方时取正号，因为被积函数恒正，积分值必然为正。这个细节看着 trivial，但在处理带符号的变形时真的要小心，后面第 7 节我会专门讲一个因为开方符号搞错导致整篇推导报废的例子。
顺带得到两个常用结论：$\int_0^{+\infty}e^{-x^2}dx=\frac{\sqrt{\pi}}{2}$，以及 $\int_{-\infty}^{+\infty}e^{-x^2/2}dx=\sqrt{2\pi}$。后者是标准正态分布归一化常数的直接来源，$e^{-x^2/2}$ 这个形式在概率论里出现频率远高于 $e^{-x^2}$，原因后面会讲。
还有一条不太常用但很漂亮的路径：用伽马函数。令 $t=x^2$，则 $dx=\frac{1}{2}t^{-1/2}dt$，于是 $\int_0^\infty e^{-x^2}dx=\frac{1}{2}\int_0^\infty t^{-1/2}e^{-t}dt=\frac{1}{2}\Gamma(\frac{1}{2})$。由 $\Gamma(\frac{1}{2})=\sqrt{\pi}$ 立刻得到同样结果。这条路的优势是可以直接推广到 $\int_0^\infty x^{n}e^{-x^2}dx$ 这类带幂次的积分——那就是 $\frac{1}{2}\Gamma(\frac{n+1}{2})$，一步到位，比反复分部积分快得多。我在算量子力学期望值的时候基本都用这个套路。
3. 大学物理里它换了四张脸，很多人没认出来是同一个东西
3.1 麦克斯韦速率分布：那个吓人的归一化常数是算出来的
麦克斯韦速率分布写成 $f(v)=4\pi\left(\frac{m}{2\pi kT}\right)^{3/2}v^2e^{-mv^2/2kT}$。前头那一坨系数是保证 $\int_0^\infty f(v)dv=1$ 的归一化因子。很多人背下来了事，其实推一遍非常痛快。
从速度分布出发，三个方向的分布是独立的，各为 $\sqrt{\frac{m}{2\pi kT}}e^{-mv_i^2/2kT}$，因为它们各自要满足一维归一化 $\int_{-\infty}^\infty e^{-mv_i^2/2kT}dv_i=\sqrt{\frac{2\pi kT}{m}}$。三个方向相乘，系数就是 $\left(\frac{m}{2\pi kT}\right)^{3/2}$，指数部分是 $e^{-m(v_x^2+v_y^2+v_z^2)/2kT}$。
再换到球坐标，$v^2=v_x^2+v_y^2+v_z^2$，速度空间体积元是 $v^2\sin\theta,dv,d\theta,d\varphi$，对角度积分给出 $4\pi$，于是得到速率分布。整个过程的本质就是：一维高斯积分负责给每个方向定系数，球坐标的雅可比因子 $v^2$ 负责把"速度向量"变成"速率标量"。
这里有个常被忽略的点：$v^2$ 这个因子不是装饰，它决定了分布形状。没有 $v^2$ 的话分布在 $v=0$ 处最大；有了 $v^2$，分布在 $v=0$ 处为零，峰值被推到某个正的速率上。这背后的物理直觉是：高速率对应的速度空间"球壳"体积更大，能装的微观状态更多，所以概率密度被"稀释"到更大的 $v$ 附近。我在给一个做分子动力学仿真的朋友解释采样问题时用过这个类比，他立刻就明白了为什么自己的初始速度采样不能直接用高斯分布而必须乘上 $v^2$。
3.2 量子谐振子基态：波函数归一化就是一维高斯积分
一维量子谐振子的基态波函数是 $\psi_0(x)=Ae^{-m\omega x^2/2\hbar}$。归一化条件 $\int_{-\infty}^\infty|\psi_0|^2dx=1$ 给出 $|A|^2\int_{-\infty}^\infty e^{-m\omega x^2/\hbar}dx=1$。套用一般公式 $\int_{-\infty}^\infty e^{-ax^2}dx=\sqrt{\pi/a}$，取 $a=m\omega/\hbar$，得到 $|A|^2\sqrt{\frac{\pi\hbar}{m\omega}}=1$，即 $A=\left(\frac{m\omega}{\pi\hbar}\right)^{1/4}$。
我在第一次自己推这个的时候，卡在了四次方根上——总以为系数应该是个平方根形式。其实仔细想就明白了：归一化条件里用的是 $|\psi|^2$，所以 $|A|^2$ 是 $\sqrt{m\omega/\pi\hbar}$，那 $A$ 自然要开一次方，得到四分之一次方。这种"归一化对象是模平方还是模本身"的混淆，我在后来做光学场强归一化时又犯过一次，教训很深。
激发态就更有意思了。$\psi_n$ 里面带着厄米多项式 $H_n$，归一化要算 $\int_{-\infty}^\infty H_m(x)H_n(x)e^{-x^2}dx=\sqrt{\pi}2^nn!\delta_{mn}$。这个正交关系本质上就是带权的高斯积分，权重函数是 $e^{-x^2}$。数值计算里对应的工具叫高斯-埃尔米特求积，后面第 6 节会详细说。量子化学里算谐振子基函数的重叠积分，天天都在用这个东西。
3.3 高斯光束与衍射包络：场强积分给出光斑尺寸
激光里的高斯光束，场强横向分布是 $E(r)=E_0e^{-r^2/w^2}$，$w$ 是束腰半径。总功率 $P=\int_0^\infty |E(r)|^2 2\pi r,dr=|E_0|^2\cdot 2\pi\int_0^\infty e^{-2r^2/w^2}r,dr$。注意这里又出现了 $r$ 这个雅可比因子，而且它和 $e^{-2r^2/w^2}$ 一配，积分直接可算：令 $u=2r^2/w^2$，得 $\int_0^\infty =\frac{w^2}{4}$，所以 $P=|E_0|^2\frac{\pi w^2}{2}$。这就是为什么说高斯光束 86.5% 的能量集中在 $r\le w$ 的圆内——这个数字来自 $\int_0^w$ 与 $\int_0^\infty$ 的比值 $1-e^{-2}\approx0.865$，背后是同一个高斯积分。
我做光学仿真的时候发现一个现象：很多人用二维数组生成高斯光斑，直接对 $e^{-(x^2+y^2)/w^2}$ 采样，然后忘了做离散归一化，导致总功率随网格分辨率变化。正确做法是要么在连续域解析归一化（乘上 $\frac{2}{\pi w^2}$），要么在离散域除以所有像素值之和。这两种做法在网格够密的时候结果一致，网格粗糙时差别明显。几年前我在一个项目里因为这个问题导致仿真出来的耦合效率总是比实验低 15%，查了两周才发现是归一化没做对。
3.4 统计物理配分函数：谐振子求和变积分时的那个 $\sqrt{\pi}$
固体物理里算声子配分函数，或者统计力学里算理想气体的平动配分函数，都会遇到 $\int_{-\infty}^\infty e^{-\beta p^2/2m}dp=\sqrt{2\pi mkT}$ 这样的式子。一维平动配分函数里动量部分的贡献就是这个积分。
更进一步，爱因斯坦模型的配分函数是 $Z=\sum_n e^{-\beta\hbar\omega(n+1/2)}$，这是求和，不是积分。但到了德拜模型、到了连续极限，求和就变成了积分，而积分核往往是高斯型的。我印象最深的是算经典极限下的谐振子配分函数：$Z=\frac{1}{h}\int\int e^{-\beta(p^2/2m+m\omega^2q^2/2)}dp,dq$，两个方向各出一个高斯积分，乘起来是 $\frac{1}{h}\sqrt{\frac{2\pi m}{\beta}}\sqrt{\frac{2\pi}{\beta m\omega^2}}=\frac{1}{h}\frac{2\pi}{\beta\omega}=\frac{kT}{\hbar\omega}$。这个结果正好是量子结果 $1/(e^{\beta\hbar\omega}-1)$ 在高温下的极限，对得上。这个验证方法我强烈推荐：算完积分后取一个已知极限检验一下，比反复检查推导过程有效得多。
4. 概率统计视角：正态分布不过是同一个积分的马甲
4.1 归一化常数 $1/\sqrt{2\pi}\sigma$ 的推导只用了三步
标准正态密度 $\varphi(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$，要验证它积分为 1，核心就是算 $\int_{-\infty}^\infty e^{-x^2/2}dx$。三步：第一，令 $x=\sqrt{2}u$，$dx=\sqrt{2}du$，积分变成 $\sqrt{2}\int e^{-u^2}du=\sqrt{2}\sqrt{\pi}=\sqrt{2\pi}$。第二，一般情形 $N(\mu,\sigma^2)$ 的密度是 $\frac{1}{\sigma\sqrt{2\pi}}e^{-(x-\mu)^2/2\sigma^2}$，令 $y=(x-\mu)/\sigma$ 就化归到标准情形。第三，验证完毕。
漂亮的地方在于，这个 $1/\sqrt{2\pi}$ 和 $\sqrt{\pi}$ 只差一个 $\sqrt{2}$，而这个 $\sqrt{2}$ 完全来自方差定义里那个习惯性的 $1/2$。如果当年有人把正态分布写成 $e^{-\pi x^2}$（事实上真有数学家这么推广过），归一化常数就变成 1，形式更简洁。之所以没这么写，纯粹是因为 $\sigma^2$ 作为方差有明确的统计意义，比"半方差"这种量更直观。
我给别人讲正态分布的时候，总会强调一点：$e^{-x^2/2}$ 里那个 $1/2$ 是方差定义带来的，$1/\sqrt{2\pi}$ 里的 $\sqrt{2\pi}$ 是高斯积分给的。这两部分来源完全独立，但凑在一起就显得很神秘。拆开看，一点都不神秘。
4.2 特征函数：把高斯积分搬到复平面上
正态分布的特征函数是 $\phi(t)=E[e^{itX}]=e^{i\mu t-\sigma^2t^2/2}$。推导过程：
$$\phi(t)=\frac{1}{\sqrt{2\pi}\sigma}\int_{-\infty}^\infty e^{itx}e^{-(x-\mu)^2/2\sigma^2}dx$$
换元 $y=(x-\mu)/\sigma$，得到 $e^{i\mu t}\cdot\frac{1}{\sqrt{2\pi}}\int e^{ity\sigma}e^{-y^2/2}dy$。剩下的积分是 $\int e^{-y^2/2+ity\sigma}dy$，这就要用带线性项的高斯积分公式（下一节会推）：$\int e^{-ay^2+by}dy=\sqrt{\pi/a},e^{b^2/4a}$，其中 $a=1/2$，$b=it\sigma$。代入得 $\sqrt{2\pi}e^{-t^2\sigma^2/2}$，再乘上前面的 $1/\sqrt{2\pi}$ 正好抵消，剩 $e^{i\mu t-\sigma^2t^2/2}$。
这条推导在概率论课本里经常被简写成"可以证明"，但它其实是高斯积分最强大的应用之一：把实轴上的积分解析延拓到复平面，让复指数和实高斯核结合。这个技巧在信号处理里就是傅里叶变换的高斯窗，在量子力学里就是相干态的相空间表示。我做过一段时间的随机过程建模，特征函数简直是神器——两个独立随机变量之和的特征函数是各自特征函数的乘积，正态分布的特征函数是 $e^{-\sigma^2t^2/2}$，两个一乘，方差直接相加，一秒钟证完正态分布的可加性，比卷积积分快了一个数量级。
4.3 多维高斯与协方差矩阵：雅可比因子换成了行列式
$d$ 维高斯分布密度是 $\frac{1}{(2\pi)^{d/2}|\Sigma|^{1/2}}\exp\left[-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x}-\boldsymbol{\mu})\right]$。要验证它积分为 1，做法是令 $\mathbf{y}=\Sigma^{-1/2}(\mathbf{x}-\boldsymbol{\mu})$，则二次型变成 $\mathbf{y}^T\mathbf{y}$，而体积元变换带出的雅可比因子是 $|\Sigma|^{1/2}$。于是积分变成 $|\Sigma|^{1/2}\int e^{-\mathbf{y}^T\mathbf{y}/2}d\mathbf{y}=|\Sigma|^{1/2}(2\pi)^{d/2}$，正好和归一化常数抵消。
这个推导里有两个点值得特别留意。第一，雅可比因子从一维的 $r$ 变成了多维的行列式，这是同一个思想的不同表现形式；第二，要求 $\Sigma$ 正定，否则 $\Sigma^{-1/2}$ 不存在，积分发散。我在做高斯过程回归的时候，协方差矩阵偶尔会因为数值问题变得接近奇异，导致 $\log|\Sigma|$ 变成负无穷，整个似然函数直接崩掉。解决办法是加一个小的对角项（jitter），也就是 $\Sigma+\epsilon I$，这个技巧在 GP 库的文档里都写着，但很少有人解释它为什么必要——本质是保证矩阵正定，让那个行列式不为零或者负。
5. 一般形式 $e^{-ax^2+bx+c}$：配方一次，全家通吃
5.1 配方三步走，别跳步
要算 $\int_{-\infty}^\infty e^{-ax^2+bx+c}dx$（要求 $a>0$），标准流程是三步。
第一步，把指数里的二次式配方：$-ax^2+bx+c=-a\left(x-\frac{b}{2a}\right)^2+\frac{b^2}{4a}+c$。这一步是纯粹的代数，但极易出错，尤其是 $b$ 是复数或者 $a$ 是矩阵的时候。
第二步，换元 $u=x-\frac{b}{2a}$，$dx=du$，积分变成 $e^{b^2/4a+c}\int e^{-au^2}du$。注意积分限不用改，因为平移不改变无穷限。
第三步，套一维高斯积分 $\int e^{-au^2}du=\sqrt{\pi/a}$（要求 $a$ 的实部为正），得到最终结果：
$$\int_{-\infty}^\infty e^{-ax^2+bx+c}dx=\sqrt{\frac{\pi}{a}}\exp\left(\frac{b^2}{4a}+c\right)$$
这个公式值得背下来，因为它是所有高斯积分变体的母公式。我个人的经验是：遇到任何高斯型积分，先想办法把它化成这个形式，再套公式，比每次重新推一遍要快得多，也不容易错。
举个实际例子：量子力学里算相干态的波函数重叠，会出 $\int e^{-\alpha x^2+\beta x}dx$ 这种，直接套用 $\sqrt{\pi/\alpha}e^{\beta^2/4\alpha}$。热统里算带外场的配分函数，也经常是这形式。我甚至在做金融里的期权定价时见过——Black-Scholes 公式里的那个正态累积分布函数，本质就是对数正态分布做高斯积分，只不过换了个变量。
5.2 一张表记住所有常见变体
我整理过一张表，把工作中真正用得上的高斯积分变体列出来，贴在这里供参考：
| 积分形式 | 结果 | 常用场景 | 备注 |
|---|---|---|---|
| $\int_{-\infty}^\infty e^{-ax^2}dx$ | $\sqrt{\pi/a}$ | 归一化、配分函数 | $a>0$，最基础形式 |
| $\int_0^\infty e^{-ax^2}dx$ | $\frac{1}{2}\sqrt{\pi/a}$ | 半无限区间 | 由对称性直接得 |
| $\int_{-\infty}^\infty e^{-ax^2+bx}dx$ | $\sqrt{\pi/a},e^{b^2/4a}$ | 特征函数、生成函数 | 配方后套用 |
| $\int_{-\infty}^\infty xe^{-ax^2}dx$ | $0$ | 奇函数期望 | 对称性，无需计算 |
| $\int_{-\infty}^\infty x^2e^{-ax^2}dx$ | $\frac{1}{2a}\sqrt{\pi/a}$ | 方差、动能期望 | 对 $a$ 求导可得 |
| $\int_0^\infty x^{2n}e^{-ax^2}dx$ | $\frac{(2n-1)!!}{2^{n+1}}\sqrt{\pi/a^{2n+1}}$ | 高阶矩、统计量 | 双阶乘形式 |
| $\int_{-\infty}^\infty e^{-ax^2}\cos(bx)dx$ | $\sqrt{\pi/a},e^{-b^2/4a}$ | 傅里叶变换 | 取实部 |
| $\int_{-\infty}^\infty e^{-ax^2}\sin(bx)dx$ | $0$ | 奇函数 | 对称性 |
这张表里最值得记的是"对参数求导"这个技巧。$\int x^2e^{-ax^2}dx$ 不需要另起炉灶，只要把 $\int e^{-ax^2}dx=\sqrt{\pi}a^{-1/2}$ 对 $a$ 求一次导，再取负号，立刻得到 $\frac{1}{2}\sqrt{\pi}a^{-3/2}$。求两次导就有 $x^4$。这个"积分号下求导"的操作在物理里特别常用，比如算能量均分定理里各自由度的贡献。不过要注意，求导和积分交换次序需要被积函数对参数一致可积，高斯型积分满足这个条件，所以可以放心操作。
5.3 带 $|x|$、带 $x^3$、带 $\sqrt{x}$ 的怪形状怎么处理
现实中的积分很少老老实实是纯高斯形式，经常冒出一些奇怪的乘子。
遇到 $|x|e^{-ax^2}$，利用偶函数性质，等于 $2\int_0^\infty xe^{-ax^2}dx=\frac{1}{a}$，非常简单，就是因为 $|x|$ 把负半轴翻到正半轴，两倍算一次。
遇到 $x^3e^{-ax^2}$，奇函数，答案是 0。这类"一眼零"的积分在物理推导里能省大量时间，前提是你养成先检查奇偶性的习惯。我见过有同学硬算了一整页纸，最后发现是零，那种心情我懂。
遇到 $\sqrt{x}e^{-ax^2}$，这个就要换元 $t=x^2$ 转成伽马函数了：$\int_0^\infty \sqrt{x}e^{-ax^2}dx=\frac{1}{2}\int_0^\infty t^{-1/4}e^{-at}dt=\frac{1}{2}a^{-3/4}\Gamma(3/4)$。伽马函数在非整数参数处没有初等表达式，只能留个 $\Gamma(3/4)$ 或者查表，或者数值算。这类积分出现在分数阶微积分、反常扩散问题的解里，工程上一般直接数值求解。
还有一类更隐蔽的：被积函数里带着阶跃函数，比如 $\int_0^\infty e^{-ax^2}dx$ 的上限是 0 到无穷而不是负无穷到正无穷。这种时候千万别习惯性套 $\sqrt{\pi/a}$，要记得取一半。我在做信号处理时对齐窗函数，因为窗口是从 0 开始而不是从负无穷开始，少除了个 2，导致频谱幅度整体偏大一倍，调试了半天才发现。
6. 数值实现：解析解拿不到的时候怎么算
6.1 高斯-埃尔米特求积：为高斯型被积函数量身定做的工具
$$\int_{-\infty}^\infty e^{-x^2}f(x)dx\approx\sum_{i=1}^n w_if(x_i)$$
其中 $x_i$ 是 $n$ 阶厄米多项式 $H_n(x)$ 的根，$w_i$ 是对应的权重。这个公式的妙处在于：对于形如 $p(x)e^{-x^2}$（$p$ 是次数不超过 $2n-1$ 的多项式）的被积函数，求积公式精确成立，没有截断误差，只有舍入误差。这意味着哪怕只用 10 个节点，算 $\int e^{-x^2}$ 就已经精确到机器精度了。
我用 Python 写过一个小脚本验证这一点，用的是numpy.polynomial.hermite_e.hermegauss，节点数取 20 的时候，算 $\int_{-\infty}^\infty e^{-x^2}dx$ 得到的结果与 $\sqrt{\pi}$ 的差值在 $10^{-15}$ 量级。相比之下，梯形法要在 $[-10,10]$ 上取几万个点才能达到 $10^{-6}$ 精度，差距不是一点半点。
不过这个工具的前提是被积函数里真的带着 $e^{-x^2}$ 这个因子。如果你的被积函数是 $e^{-ax^2}$ 或者 $e^{-(x-\mu)^2/2\sigma^2}$，得先做变量替换化归到标准形式，替换过程中要注意雅可比因子别漏。具体来说，$\int e^{-(x-\mu)^2/2\sigma^2}f(x)dx$ 换元 $t=(x-\mu)/(\sqrt{2}\sigma)$，变成 $\sqrt{2}\sigma\int e^{-t^2}f(\sqrt{2}\sigma t+\mu)dt$，然后才能套高斯-埃尔米特。我有一次忘了乘 $\sqrt{2}\sigma$，结果差了将近 3 倍，排查了一晚上。
6.2 截断上限怎么定：从 $10^{-16}$ 倒推
如果不想用高斯-埃尔米特，就用最笨也最稳的办法：把积分区间截断到 $[-L,L]$，然后用辛普森法或者scipy.integrate.quad。问题是 $L$ 取多大合适。
我的经验公式是：$e^{-aL^2}\le\epsilon$，取自然对数，得到 $L\ge\sqrt{-\ln\epsilon/a}$。对于 $\epsilon=10^{-15}$（双精度机器精度）、$a=1$，得到 $L\ge\sqrt{15\ln10}\approx\sqrt{34.5}\approx5.87$。取 $L=6$ 就足够了。对于 $e^{-x^2/2}$（$a=1/2$），$L\ge\sqrt{2\times34.5}\approx8.3$，取 $L=9$。
这个估算在实际代码里非常好用。我见过不少人图省事直接取 $L=100$ 或者 $L=1000$，结果积分点数暴增，速度慢十倍，精度一点没提高，反而因为点数太多引入更多舍入误差。原理就在这：被积函数在 $|x|>L$ 之后已经小于机器精度，你再取更远的地方只是浪费算力。
另一个常见问题是对于带参数的积分，比如 $\int e^{-ax^2}dx$，$a$ 很小的时候（比如 $a=0.01$）需要 $L\approx58$，$a$ 很大的时候只需要很小的区间。写成自适应形式比较稳妥：根据 $a$ 动态计算 $L$。我通常会用a的取值范围先做一次批量估算，避免每次调用都重新算。
6.3 一份可以直接抄的 Python 验证代码
import numpy as np from scipy.integrate import quad from scipy.special import erf import numpy.polynomial.hermite_e as H # 基础验证：int exp(-x^2) dx = sqrt(pi) exact = np.sqrt(np.pi) # 方法1：自适应数值积分 v1, err1 = quad(lambda x: np.exp(-x**2), -np.inf, np.inf) # 方法2：高斯-埃尔米特求积（20节点） nodes, weights = H.hermegauss(20) v2 = np.sum(weights * np.exp(-nodes**2) * np.exp(nodes**2)) # 注意权重已含e^{-x^2} # 更标准的写法是直接把 f(x)=1 代入 v2 = np.sum(weights) # 方法3：截断 + 辛普森法 L = 6.0 xs = np.linspace(-L, L, 20001) v3 = np.trapz(np.exp(-xs**2), xs) print(f"exact = {exact:.15f}") print(f"quad = {v1:.15f} err={err1:.2e}") print(f"hermegauss= {v2:.15f}") print(f"trapz L=6 = {v3:.15f}") # 一般形式验证：int exp(-a x^2 + b x + c) dx a, b, c = 2.5, 1.3, -0.4 exact_g = np.sqrt(np.pi/a) * np.exp(b**2/(4*a) + c) v4, _ = quad(lambda x: np.exp(-a*x**2 + b*x + c), -np.inf, np.inf) print(f"general exact = {exact_g:.12f}, quad = {v4:.12f}") # 误差函数验证：int_0^x exp(-t^2) dt = sqrt(pi)/2 * erf(x) x0 = 1.7 left, _ = quad(lambda t: np.exp(-t**2), 0, x0) right = np.sqrt(np.pi)/2 * erf(x0) print(f"erf check: {left:.12f} vs {right:.12f}")这段代码我反复用过，有几个地方值得说。hermegauss返回的权重已经包含了 $e^{-x^2}$，所以你只需要把 $f(x)$ 代进去求和，别再乘一次 $e^{-x^2}$，这是最常见的误用。另外quad直接支持-np.inf到np.inf，内部会做变量替换到有限区间，不用自己截断。
还有一个细节：trapz在 L=6、两万点的时候，精度大约是 $10^{-11}$ 量级，够用但不是最好。如果精度要求高，可以把 L 加大或者换成quad。我在实际项目里的选择原则是：被积函数光滑、区间有限，优先quad；需要反复调用、追求速度，就用固定节点的高斯-埃尔米特或者预先算好的查找表。
7. 踩过的坑：这五个错误我保证你至少会犯一个
7.1 交换积分次序想当然
前面提过，$I^2$ 拆成二重积分需要绝对收敛。我第一次自己推的时候完全没意识到这是个条件，直接就拆了。后来做傅里叶变换相关的推导，遇到条件收敛的积分，照葫芦画瓢地拆，结果得出一个自相矛盾的结论，查了半天才想起来富比尼定理的存在。
判断方法其实不难：看绝对值积分是否有限。对于高斯型，$|e^{-(x^2+y^2)}|=e^{-(x^2+y^2)}$，积分有限，放心拆。对于带振荡因子的，比如 $e^{-ax^2}\sin(bx)$，绝对值积分 $\int|e^{-ax^2}\sin(bx)|dx$ 也是有界的，因为 $e^{-ax^2}$ 衰减够快。真正危险的是像 $\sin x/x$ 那种衰减慢的。我的建议是：凡是要交换次序或者拆分积分，先花十秒钟检查一下绝对可积性，这个习惯能救你无数次。
7.2 极坐标换元漏掉雅可比因子 $r$
这个前面已经强调过，但值得再提一次，因为它太容易漏了。在手工推导时，写着写着就把 $dxdy$ 写成 $drd\theta$，潜意识里觉得"换个坐标而已，面积元应该不变"。实际上除了极坐标，球坐标也有雅可比 $r^2\sin\theta$，柱坐标有 $r$，各种曲线坐标都有自己的雅可比。
我的记忆方法是：极坐标看"扇环"，球坐标看"楔形"。极坐标的扇环面积约等于 $r\cdot dr\cdot d\theta$，所以雅可比是 $r$；球坐标的楔形体积约等于 $(r\sin\theta,d\varphi)\cdot(rd\theta)\cdot dr=r^2\sin\theta,dr,d\theta,d\varphi$，所以雅可比是 $r^2\sin\theta$。用几何图像记，比死记公式靠谱。
7.3 开方符号处理不当
$\sqrt{I^2}=|I|$ 而不是 $I$，这是初中知识，但我在推导中还是栽过。有一次算一个带符号的积分，中间步骤平方了，最后开方时想当然地取了正号，结果和物理直觉矛盾。回头看才发现，被积函数在某个区间是负的，积分结果应该是负的，我取了正号就等于把符号信息丢了。
防范方法很简单：每次从 $I^2$ 回到 $I$，都明确写下 $I=\pm\sqrt{I^2}$，然后用物理或者几何意义确定符号。对于高斯积分本身，被积函数恒正，取正号没问题，但如果你在处理带相位的复积分，或者被积函数有零点穿越，就必须小心。
7.4 形状像高斯就套公式
这是最隐蔽的一个坑。比如 $\int_{-\infty}^\infty e^{-|x|}dx$，形状也像个钟形，但它不是高斯形式，套 $\sqrt{\pi/a}$ 会得到完全错误的结果（正确答案是 2）。再比如 $\int e^{-x^4}dx$，这也衰减得很快，但它不是高斯积分，没有 $\sqrt{\pi/a}$ 这种简单形式，只能数值算或者用伽马函数（$\int_0^\infty e^{-x^4}dx=\Gamma(5/4)$）。
区分方法就是看指数是不是二次型。是二次的，才能套高斯公式。是 $|x|$、$x^4$、$x^{2/3}$，都不行。我在做信号处理时见过有人把拉普拉斯分布的密度 $e^{-|x|}$ 当成高斯来处理，算出来的方差明显偏小，模型选择和实际情况差很远。
7.5 数值实现里忘了离散归一化
最后这个坑比较工程化。用数值方法生成高斯函数时，连续域的归一化常数 $\frac{1}{\sqrt{2\pi}\sigma}$ 在离散网格上不一定保证求和为 1。原因是梯形法或者中点法的离散近似和连续积分有偏差，网格越粗偏差越大。
正确的做法是：生成离散样本后，显式除以样本和（或者样本和乘以网格间距，取决于你用的是密度还是概率质量）。我在一个做贝叶斯推断的项目里，因为忘了这一步，导致后验分布的总概率不等于 1，算出来的置信区间全都偏窄，模型评估指标自然一塌糊涂。修复方式就是在采样之后加一行p = p / p.sum()（如果采样点在网格上且等间距，这行等价于除以网格间距乘以和），一秒钟的事，但没意识到的时候能卡好几天。
最后说点个人体会
高斯积分在我眼里是数学里少有的"一鱼多吃"的典型。它把高等数学的二重积分技巧、大学物理的归一化需求、概率统计的分布理论串在了一起，而这三者用的是同一个核心结论：$\int e^{-x^2}dx=\sqrt{\pi}$。真正理解它，不是记住这个公式，而是理解为什么会有 $\sqrt{\pi}$、什么时候能套、什么时候不能套、数值实现时误差从哪里来。
我自己学它的路径是：先在高数课上背结论，然后在物理课上反复看到它却认不出来，最后在写代码时因为参数错误被迫从头推一遍。推完之后，之前那些零散的公式突然就串成一条线了。如果你现在正卡在某个归一化常数上，我的建议是别急着背，找个安静的晚上自己从 $I^2$ 开始推一遍，推完你会发现自己对二重积分、极坐标、雅可比因子的理解都会上一个台阶。
还有一个小技巧分享：遇到任何高斯型积分，先在纸上写下母公式 $\int e^{-ax^2+bx+c}dx=\sqrt{\pi/a},e^{b^2/4a+c}$，然后把你的积分往这个形式凑，凑不出来就说明它不是高斯型，老老实实上数值。这个检查步骤花不了三十秒，但能帮你省掉后面一大堆返工。
