5.5 一元积分学与微积分基本定理

微分回答的是「局部怎么变」，积分回答的是「整体累积了多少」。把一条曲线下的面积切成一列窄矩形再求和，这是积分最初的想法；把总位移看成速度对时间的累积，这是积分的物理含义。上一节末尾出现的积分余项 $R_n(x)={\int }_{x_0}^x\frac{f^{(n+1)}(t)}{n!}(x-t{)}^{n}dt$ 也说明，积分的语言在证明里同样不可缺少。

本节的组织沿三条线索展开。定义线索是 Riemann 积分：分割、求和的极限，可积性的判别（Darboux 上下和），以及可积函数类的刻画。定理线索是微积分基本定理：变上限积分的导数等于被积函数，由此得到 Newton-Leibniz 公式，把求积分化为求原函数，微分与积分从此互为逆运算。推广线索把积分推向两类反常情形：无穷区间与无界函数，判敛工具是比较判别法，条件收敛的例子来自 ${\int }_1^{\mathrm{\infty }}\frac{\sin x}{x}dx$，Gamma 函数则同时带有两类反常。最后回到计算：原函数可用初等函数表示的情形在应用中是少数，数值求积承担了大部分实际计算，本节实测中点、梯形、Simpson 与 Gauss 求积的误差阶。

5.5.1 Riemann 积分

分割、Riemann 和与定义

设 $f$ 是 $[a,b]$ 上的有界函数。取分割

$$P:a=x_0<x_1<\cdots <x_n=b$$

记小区间长度 $\mathrm{\Delta }{x}_i=x_i-x_{i-1}$，在每个小区间上任取一点 ${\xi }_i\in [x_{i-1},x_i]$，作 Riemann 和

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

令分割的细度 $\lambda (P)=\underset{i}{max}\mathrm{\Delta }{x}_i$。若当 $\lambda (P)\to 0$ 时 $S(P,\xi )$ 趋于一个与分割方式和取点方式都无关的有限数 $I$，称 $f$ 在 $[a,b]$ 上 Riemann 可积，$I$ 是 $f$ 在 $[a,b]$ 上的定积分：

$$I={\int }_a^{b}f(x)dx=\underset{\lambda (P)\to 0}{lim}\sum _{i=1}^{n}f({\xi }_i)\mathrm{\Delta }{x}_i$$

Riemann 和的几何含义是「一列窄矩形面积之和」，积分的几何含义是曲边梯形的面积。当 $f$ 取负值时面积带负号，积分是带符号的累积量；物理含义是累积：速度对时间积分得位移，功率对时间积分得能量。

Darboux 上下和与可积判别

定义中的极限要求对所有取点方式成立，直接验证很麻烦。Darboux 用最值的上下界把取点的自由度消掉：记 $M_i=\underset{[x_{i-1},x_i]}{sup}f$，$m_i=\underset{[x_{i-1},x_i]}{inf}f$，作上和与下和

$$U(P)=\sum _{i=1}^n{M}_i\mathrm{\Delta }{x}_i,L(P)=\sum _{i=1}^n{m}_i\mathrm{\Delta }{x}_i$$

任何 Riemann 和都夹在 $L(P)$ 与 $U(P)$ 之间。可积的判别准则（Darboux 准则）：$f$ 在 $[a,b]$ 上有界且可积，等价于

$$\mathrm{\forall }\epsilon >0,\mathrm{\exists }P,U(P)-L(P)<\epsilon$$

即上下和之差可以任意小。证明的要点是加细分割时上和不增、下和不减，两者的极限（上积分与下积分）存在，可积等价于两个极限相等。这条准则把「极限与取点无关」的要求替换成「一种分割就能把振幅压小」的可操作条件，后续所有可积性的证明都从它出发。

可积函数类

由 Darboux 准则可以证明三类函数的可积性：

1. 连续函数：闭区间上连续则一致连续（5.3 节），把区间取得足够细可使每个小区间上的振幅 $M_i-m_i$ 一致地小于 $\frac{\epsilon }{b-a}$，于是 $U-L<\epsilon$。
2. 单调函数：单调函数在每个小区间上的振幅总和被总变差 $|f(b)-f(a)|$ 控制，$U-L\le |f(b)-f(a)|\underset{i}{max}\mathrm{\Delta }{x}_i$，细度足够小即可。
3. 有界且间断点有限的函数：把间断点附近的小区间单独拿出来（总长度可以任意小），其余部分按连续函数处理。

更一般的刻画是 Riemann-Lebesgue 定理：有界函数 Riemann 可积，当且仅当它的间断点集是零测集（零测集指可以用总长度任意小的区间覆盖的集合，5.10 节会给出严格定义）。这个定理说明间断点可以有无穷多个，只要它们足够稀少。相应的反面例子是 Dirichlet 函数：有理点取 $1$、无理点取 $0$。任何区间内都同时含有有理数与无理数，因此每个小区间上 $M_i=1$、$m_i=0$，$U(P)-L(P)=b-a$ 对任何分割都等于区间长度，不可积。

三条常用性质：改变有限个点的函数值不影响积分值（零测集上的差异不改变上下和）；$|f|$ 可积但反过来不成立；$f$ 可积则 $|f|$ 可积，且 $\left|{\int }_a^{b}f\right|\le {\int }_a^b|f|$。

5.5.2 微积分基本定理

变上限积分

设 $f$ 在 $[a,b]$ 上可积，定义变上限积分

$$F(x)={\int }_a^{x}f(t)dt,x\in [a,b]$$

两条性质把积分与微分接通。

连续性：$f$ 有界（$|f|\le M$），则对 $x,y\in [a,b]$，

$$|F(x)-F(y)|=|{\int }_y^{x}f(t)dt|\le M|x-y|$$

$F$ 满足 Lipschitz 条件，从而连续。这条估计的写法在后面反复出现：把函数的差表示成积分，再用被积函数的界控制。

可导性：若 $f$ 在点 $x_0$ 连续，则 $F$ 在 $x_0$ 可导，且

$$F^{\prime }(x_0)=f(x_0)$$

证明用差商的估值：$\frac{F(x_0+h)-F(x_0)}{h}-f(x_0)=\frac{1}{h}{\int }_{x_0}^{x_0+h}[f(t)-f(x_0)]dt$，由 $f$ 在 $x_0$ 连续，被积函数在 $t$ 靠近 $x_0$ 时一致地小，右端趋于零。结论有两个读法：一是连续函数必有原函数（$F$ 就是一个原函数），二是求导与积分在连续函数类上互为逆运算。

与微分学中值定理的对照

5.4 节的积分余项 $R_n(x)={\int }_{x_0}^x\frac{f^{(n+1)}(t)}{n!}(x-t{)}^{n}dt$ 正是「把余项写成积分再用估值控制」的用法；本节的可导性证明与它有相同的结构：差商写成积分，被积函数的振幅小则结果小。这条技巧在数值分析、偏微分方程与最优化中都以不同形式出现。

Newton-Leibniz 公式

设 $f$ 在 $[a,b]$ 上连续，$G$ 是 $f$ 的任一原函数（$G^{\prime }=f$），则

$${\int }_a^{b}f(x)dx=G(b)-G(a)\text{defeq}G(x){|}_a^b$$

这就是微积分基本定理（Newton-Leibniz 公式）。证明由「原函数相差常数」得到：$F$ 与 $G$ 都是 $f$ 的原函数，$G-F$ 的导数为零，因而是常数 $C$；取 $x=a$ 得 $F(a)=0$、$C=G(a)$；取 $x=b$ 得 ${\int }_a^{b}f=F(b)=G(b)-G(a)$。

公式的意义在于计算方式的转换：积分的定义是「分割求和的极限」，直接计算需要处理无穷过程；有了原函数，一次求值相减即可得到结果。代价是求原函数并不总有初等表达：$\frac{\sin x}{x}$、$e^{-x^2}$、$\frac{1}{\ln x}$ 的原函数都超出初等函数的范围（可以证明），此时积分依然存在（被积函数连续），但只能数值计算，这正是 5.5.5 与实验一处理的情形。

两个基本技巧是求原函数的常规工具。

换元积分：

$${\int }_a^{b}f(g(t))g^{\prime }(t)dt={\int }_{g(a)}^{g(b)}f(u)du$$

方向是「把复杂的被积函数通过替换化简」。例：${\int }_0^1\frac{2x}{1+x^2}dx$，取 $u=1+x^2$、$du=2xdx$，原式 $={\int }_1^2\frac{du}{u}=\ln 2=0.693147\dots$。

分部积分：

$${\int }_a^{b}udv=uv{|}_a^b-{\int }_a^{b}vdu$$

来源是乘积求导法则 $(uv{)}^{\prime }=u^{\prime }v+u{v}^{\prime }$ 的移项积分。例：${\int }_0^{1}x{e}^{-x}dx$，取 $u=x$、$dv=e^{-x}dx$，得

$${\int }_0^{1}x{e}^{-x}dx=-x{e}^{-x}{|}_0^1+{\int }_0^1{e}^{-x}dx=-e^{-1}+(1-e^{-1})=1-\frac{2}{e}=0.2642411177\dots$$

这个值在练习题中还会以数值求积的方式复算一次。分部积分在理论上的用途更广：Gamma 函数的递推公式、含参积分的估计、弱导数的定义都由它导出。

5.5.3 反常积分

Riemann 积分要求区间有界、函数有界。两类越界的情形可以分别定义极限来救回：无穷区间与无界函数。

无穷区间上的积分

设 $f$ 在 $[a,+\mathrm{\infty })$ 的每个有限子区间上可积，定义

$${\int }_a^{+\mathrm{\infty }}f(x)dx=\underset{R\to +\mathrm{\infty }}{lim}{\int }_a^{R}f(x)dx$$

右端极限存在（有限）时称积分收敛，否则称发散。典型例子是 $p$ 积分：

$${\int }_1^{+\mathrm{\infty }}\frac{dx}{x^p}\text{收敛}⟺p>1,{\int }_1^{+\mathrm{\infty }}\frac{dx}{x^p}=\frac{1}{p-1}(p>1)$$

$p\le 1$ 时发散：$p=1$ 的部分积分是 $\ln R$（对数发散），$p<1$ 的部分积分按 $R^{1-p}$ 增长（幂律发散）。判据的含义是衰减速度：被积函数趋于零还不够，必须快于 $\frac{1}{x}$。

无界函数的积分

设 $f$ 在 $b$ 的左侧无界（奇点），定义

$${\int }_a^{b}f(x)dx=\underset{\epsilon \to 0^{+}}{lim}{\int }_a^{b-\epsilon }f(x)dx$$

典型例子是 $p$ 积分在零点附近的行为：

$${\int }_0^1\frac{dx}{x^p}\text{收敛}⟺p<1$$

阈值与无穷区间的情形恰好相反。原因在于两类反常的尺度不同：无穷远处考察的是衰减是否快于 $\frac{1}{x}$（$p>1$ 才够快），零点附近考察的是增长是否慢于 $\frac{1}{x}$（$p<1$ 才够慢）。同一个函数 $\frac{1}{x}$ 在两个位置都处于临界状态：${\int }_1^{\mathrm{\infty }}\frac{dx}{x}$ 与 ${\int }_0^1\frac{dx}{x}$ 都发散。

比较判别法与条件收敛

判断收敛不需要算出积分值，用已知的 $p$ 积分作比较即可。比较判别法：设 $0\le f(x)\le g(x)$ 在 $[a,+\mathrm{\infty })$ 上成立，若 ${\int }_a^{\mathrm{\infty }}g$ 收敛则 ${\int }_a^{\mathrm{\infty }}f$ 收敛（充要配对：$f\ge g$ 且 $\int g$ 发散则 $\int f$ 发散）。极限形式：若 $\frac{f(x)}{g(x)}\to L$（$0<L<+\mathrm{\infty }$），则两个积分同敛散。例：${\int }_1^{\mathrm{\infty }}\frac{x+1}{x^3+2}dx$，被积函数与 $\frac{1}{x^2}$ 同阶（比值趋于 $1$），$p=2>1$，收敛。

绝对收敛与条件收敛：$\int |f|$ 收敛称绝对收敛，绝对收敛蕴含收敛（用 $|f|$ 的收敛控制 $f$ 的部分积分满足 Cauchy 准则）。收敛但不绝对收敛称条件收敛，标准例子是

$${\int }_1^{+\mathrm{\infty }}\frac{\sin x}{x}dx\text{条件收敛}$$

它的收敛由 Dirichlet 判别法保证：被积函数写成 $a(x)b(x)$，$a$ 的积分族有界（$|{\int }_1^R\sin xdx|\le 2$），$b(x)=\frac{1}{x}$ 单调趋于零。两个条件分别提供「部分积分不累积」与「振荡幅度衰减」，合起来保证部分积分满足 Cauchy 准则；绝对发散是因为 $|\sin x|$ 的平均值约为 $\frac{2}{\pi }$，$\int \frac{|\sin x|}{x}dx$ 与 $\int \frac{dx}{x}$ 同阶。

一个提醒：Cauchy 主值与收敛是两件事。对称截断 $\underset{R\to \mathrm{\infty }}{lim}{\int }_{-R}^{R}xdx=0$ 存在，但 ${\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}xdx$ 不收敛（两侧分别发散）。条件收敛的积分对截断方式敏感，数值计算时要按定义双侧分别取极限。

5.5.4 Gamma 与 Beta 函数

Gamma 函数把阶乘推广到连续参数：

$$\mathrm{\Gamma }(s)={\int }_0^{+\mathrm{\infty }}{t}^{s-1}{e}^{-t}dt,s>0$$

这个积分同时带有两类反常：$t\to 0^{+}$ 处 $t^{s-1}$ 的无界性要求 $s>0$（用 $p$ 积分判据，$1-s<1$），$t\to +\mathrm{\infty }$ 处指数衰减保证收敛（任意指数衰减都快于任何幂函数）。

两条基本性质：

递推：由分部积分（取 $u=t^s$、$dv=e^{-t}dt$），

$$\mathrm{\Gamma }(s+1)={\int }_0^{\mathrm{\infty }}{t}^s{e}^{-t}dt=-t^s{e}^{-t}{|}_0^{\mathrm{\infty }}+s{\int }_0^{\mathrm{\infty }}{t}^{s-1}{e}^{-t}dt=s\mathrm{\Gamma }(s)$$

结合 $\mathrm{\Gamma }(1)=1$ 得 $\mathrm{\Gamma }(n+1)=n!$，Gamma 函数是阶乘的唯一对数凸延拓（Bohr-Mollerup 定理）。

特殊值 $\mathrm{\Gamma }\left(\frac{1}{2}\right)=\sqrt{\pi }$。证明用换元 $t=u^2$：${\int }_0^{\mathrm{\infty }}{t}^{-1/2}{e}^{-t}dt=2{\int }_0^{\mathrm{\infty }}{e}^{-u^2}du$，后者是 Gauss 积分的一半，值为 $\sqrt{\pi }$。这条换元在实验三中会再算一次。

Beta 函数定义为

$$B(p,q)={\int }_0^1{t}^{p-1}(1-t{)}^{q-1}dt,p,q>0$$

它与 Gamma 函数的关系是 $B(p,q)=\frac{\mathrm{\Gamma }(p)\mathrm{\Gamma }(q)}{\mathrm{\Gamma }(p+q)}$（证明用换元把两类积分拼接）。三者构成组合计算的基本工具箱：$\mathrm{\Gamma }(n+1)=n!$ 处理排列，$B$ 处理二项式型的积分，$\frac{1}{B(p,q)}$ 出现在 Beta 分布中。

Stirling 公式给出阶乘的渐近展开：

$$n!\sim \sqrt{2\pi n}{\left(\frac{n}{e}\right)}^n,\ln n!=n\ln n-n+\frac{1}{2}\ln (2\pi n)+\frac{1}{12n}-\frac{1}{360n^3}+O(n^{-5})$$

第一项之后的修正是 $\frac{1}{12n}$，实验二的数据会显示「$n$ 乘以误差」稳定在 $\frac{1}{12}=0.0833$。Stirling 公式的应用范围很广：信息论中推导熵与组合数的指数率（第四章），算法分析中估计阶乘规模的数据结构代价（例如排列枚举、哈希分析），统计学中近似二项分布与多项分布的系数（第二章）。

与概率论一章的连接

${\int }_0^{\mathrm{\infty }}{e}^{-t}{t}^{s-1}dt$ 的形状正是 Gamma 分布的归一化常数：密度 $\frac{1}{\mathrm{\Gamma }(k){\theta }^k}{x}^{k-1}{e}^{-x/\theta }$ 的积分等于 1 直接来自 Gamma 函数的定义。指数分布、${\chi }^2$ 分布、Beta 分布都是它的变体，2.5 节的连续型分布部分会用到本节全部结论。

5.5.5 数值求积

原函数可用初等函数表示的情形在应用中占少数，数值求积（quadrature）承担其余的计算。基本做法是用简单的近似积分替换原积分：在区间上取若干节点，用节点值的加权和作为积分近似。

插值型求积：用 $f$ 在节点上的插值多项式代替 $f$ 再积分。三种常用格式（单区间 $[a,b]$，$m=\frac{a+b}{2}$）：

中点公式：${\int }_a^{b}f\approx (b-a)f(m)$，误差 $=\frac{(b-a{)}^3}{24}{f}^{″}(\xi )$。

梯形公式：${\int }_a^{b}f\approx \frac{b-a}{2}[f(a)+f(b)]$，误差 $=-\frac{(b-a{)}^3}{12}{f}^{″}(\xi )$。

Simpson 公式：${\int }_a^{b}f\approx \frac{b-a}{6}[f(a)+4f(m)+f(b)]$，误差 $=-\frac{(b-a{)}^5}{2880}{f}^{(4)}(\xi )$。

中点与梯形的误差系数是 $2$ 倍关系（$\frac{1}{24}$ 与 $\frac{1}{12}$），但符号相反，这个「反号」使两者可以组合出精确度更高的公式。Simpson 公式用三个节点达到三阶代数精度（对三次多项式精确），误差阶由 $f^{(4)}$ 控制。

复化格式：把 $[a,b]$ 等分成 $n$ 段（步长 $h=\frac{b-a}{n}$），在每个小段上用上面的公式再求和。复化误差是

$$\text{中点}:\frac{(b-a)h^2}{24}{f}^{″}(\xi ),\text{梯形}:-\frac{(b-a)h^2}{12}{f}^{″}(\xi ),\text{Simpson}:-\frac{(b-a)h^4}{180}{f}^{(4)}(\xi )$$

误差阶分别是 $O(h^2)$、$O(h^2)$、$O(h^4)$：段数加倍，中点与梯形的误差降到四分之一，Simpson 降到十六分之一。这是实验一实测的入口。

Gauss 求积走另一条路线：节点不预先等距，而是与权重一起挑选，让误差的自由度最大。$n$ 个节点的 Gauss-Legendre 公式

$${\int }_{-1}^{1}f(x)dx\approx \sum _{i=1}^n{w}_{i}f(x_i)$$

代数精度达到 $2n-1$（能精确积分所有次数不超过 $2n-1$ 的多项式，是 $n$ 节点公式的理论上限）。节点是 Legendre 多项式 $P_n$ 的零点，两点公式是一个常用的手算例子：

$${\int }_{-1}^{1}f(x)dx\approx f(-\frac{1}{\sqrt{3}})+f\left(\frac{1}{\sqrt{3}}\right)$$

对光滑函数，Gauss 求积的误差随 $n$ 增加以超几何速度下降：实验一中 $n=5$ 个节点把 $[0,1]$ 上 $e^x$ 的积分误差压到 $6.5\times {10}^{-13}$，比复化 Simpson 用 128 段（129 次求值）还低两个数量级，而只用了 5 次求值。代价是节点与权重需要预计算（一般用特征值方法或递推生成），以及误差分析依赖被积函数的解析性。含端点奇性的积分（实验三的 $x^{-1/2}$）不满足解析性假设，需要先换元或用 Gauss-Jacobi 型公式。

选型经验：被积函数光滑、需要高精度、求值代价高（一次函数求值对应一次仿真或一次实验测量）时选 Gauss 求积；函数在网格上采样得到、或含有轻微的端点效应时用复化 Simpson（无奇性）与复化梯形（周期函数，谱精度）；只需要低精度快速估计时用中点或梯形，并配合 Richardson 外推提高两阶精度。

5.5.6 数值实验

实验一：求积公式的误差阶

计算 $I={\int }_0^1{e}^{x}dx=e-1$，观察四种公式的误差随节点数的变化。三种复化公式的段数取 $n=2,4,\dots ,128$，「阶」列是相邻两次误差比值的以 2 为底的对数，对应「段数加倍时误差下降的幂次」。

import math
import numpy as np

def midpoint(f, a, b, n):
    x = a + (np.arange(n) + 0.5) * (b - a) / n
    return (b - a) / n * float(np.sum(f(x)))

def trapezoid(f, a, b, n):
    x = np.linspace(a, b, n + 1)
    y = f(x)
    return (b - a) / n * float(0.5 * (y[0] + y[-1]) + np.sum(y[1:-1]))

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

def gauss(f, a, b, n):
    x, w = np.polynomial.legendre.leggauss(n)
    t = 0.5 * (b - a) * x + 0.5 * (a + b)
    return 0.5 * (b - a) * float(np.sum(w * f(t)))

f = np.exp
exact = math.e - 1.0
print(f"I = int_0^1 e^x dx = e - 1 = {exact:.15f}")
ns = [2, 4, 8, 16, 32, 64, 128]
err_mid, err_trp, err_smp = [], [], []
print(f"{'n':>5} {'中点误差':>13} {'阶':>5} {'梯形误差':>13} {'阶':>5} {'Simpson 误差':>13} {'阶':>5}")
for i, n in enumerate(ns):
    em = abs(midpoint(f, 0, 1, n) - exact)
    et = abs(trapezoid(f, 0, 1, n) - exact)
    es = abs(simpson(f, 0, 1, n) - exact)
    err_mid.append(em)
    err_trp.append(et)
    err_smp.append(es)
    om = f"{math.log2(err_mid[i-1] / em):.2f}" if i else "-"
    ot = f"{math.log2(err_trp[i-1] / et):.2f}" if i else "-"
    os_ = f"{math.log2(err_smp[i-1] / es):.2f}" if i else "-"
    print(f"{n:>5} {em:>13.4e} {om:>5} {et:>13.4e} {ot:>5} {es:>13.4e} {os_:>5}")

print("Gauss-Legendre 求积（n 个节点，不分段）")
for n in range(1, 9):
    e = abs(gauss(f, 0, 1, n) - exact)
    print(f"{n:>5} {e:>13.4e}")

输出：

I = int_0^1 e^x dx = e - 1 = 1.718281828459045
    n          中点误差     阶          梯形误差     阶    Simpson 误差     阶
    2    1.7769e-02     -    3.5649e-02     -    5.7932e-04     -
    4    4.4665e-03  1.99    8.9401e-03  2.00    3.7013e-05  3.97
    8    1.1182e-03  2.00    2.2368e-03  2.00    2.3262e-06  3.99
   16    2.7964e-04  2.00    5.5930e-04  2.00    1.4559e-07  4.00
   32    6.9915e-05  2.00    1.3983e-04  2.00    9.1027e-09  4.00
   64    1.7479e-05  2.00    3.4958e-05  2.00    5.6897e-10  4.00
  128    4.3698e-06  2.00    8.7396e-06  2.00    3.5562e-11  4.00
Gauss-Legendre 求积（n 个节点，不分段）
    n            误差
    1    6.9561e-02
    2    3.8545e-04
    3    8.2409e-07
    4    9.3297e-10
    5    6.5370e-13
    6    0.0000e+00
    7    2.2204e-16
    8    2.2204e-16

三段数据分别对应三种公式的误差阶。中点与梯形的阶列稳定在 $2.00$，与 $O(h^2)$ 一致；误差的系数也与理论吻合：$n=128$ 时 $h=\frac{1}{128}$，中点误差的理论量级 $\frac{h^2}{24}e\approx 6.9\times {10}^{-6}$，实测 $4.37\times {10}^{-6}$（系数中的 $f^{″}(\xi )=e^{\xi }$ 取区间内某点，比 $e$ 小）。中点误差与梯形误差的比值保持在 $2$ 附近（$1.7769\times {10}^{-2}$ 与 $3.5649\times {10}^{-2}$ 之比是 $2.01$），接近理论系数比 $\frac{1/24}{1/12}$ 的 $2$ 倍，符号相反在数值上表现为「中点偏小、梯形偏大」交替出现。

Simpson 的阶列从 $3.97$ 稳定到 $4.00$，与 $O(h^4)$ 一致，误差从 $n=2$ 的 $5.79\times {10}^{-4}$ 降到 $n=128$ 的 $3.56\times {10}^{-11}$，六次加倍共降低 $\frac{5.79\times {10}^{-4}}{3.56\times {10}^{-11}}\approx 1.6\times {10}^7\approx 4^6$（每加倍一次误差降到十六分之一）。

Gauss 求积的误差下降速度明显快于任何固定阶：$n=1$ 时 $6.96\times {10}^{-2}$，$n=4$ 时 $9.33\times {10}^{-10}$，$n=5$ 时 $6.54\times {10}^{-13}$。$n\ge 6$ 之后误差降到机器精度以下，表中出现精确的零（$n=6$ 的结果与 $e-1$ 的双精度值完全一致）与 $2.22\times {10}^{-16}$（一个机器精度单位）这样的舍入噪声，不再有阶的含义。作为对照：达到 ${10}^{-13}$ 量级，复化 Simpson 需要约 ${10}^5$ 次求值（按 $h^4$ 外推），Gauss 需要 5 次，区别在于 Gauss 的节点分布（向端点加密、内部稀疏）恰好匹配光滑函数的误差分布。

实验二：Stirling 近似的精度阶梯

比较 $\ln (n!)$ 的精确值与 Stirling 近似，观察误差如何随 $n$ 变化。精确值用 math.lgamma（Gamma 函数的对数，内部实现避免了 $n!$ 溢出的问题）。

import math

print(f"{'n':>6} {'ln(n!) 精确值':>16} {'Stirling':>16} {'误差':>11} {'n*误差':>10} "
      f"{'加 1/(12n) 修正':>16} {'误差':>11}")
for n in [5, 10, 20, 50, 100, 200, 500, 1000]:
    ex = math.lgamma(n + 1)
    st = n * math.log(n) - n + 0.5 * math.log(2 * math.pi * n)
    st2 = st + 1.0 / (12.0 * n)
    print(f"{n:>6} {ex:>16.12f} {st:>16.12f} {abs(st - ex):>11.3e} {n * abs(st - ex):>10.4f} "
          f"{st2:>16.12f} {abs(st2 - ex):>11.3e}")

输出：

     n       ln(n!) 精确值         Stirling          误差       n*误差     加 1/(12n) 修正          误差
     5   4.787491742782   4.770847051592   1.664e-02     0.0832   4.787513718259   2.198e-05
    10  15.104412573076  15.096082009642   8.331e-03     0.0833  15.104415342975   2.770e-06
    20  42.335616460753  42.331450141061   4.166e-03     0.0833  42.335616807728   3.470e-07
    50 148.477766951773 148.476100307326   1.667e-03     0.0833 148.477766973993   2.222e-08
   100 363.739375555563 363.738542225008   8.333e-04     0.0833 363.739375558341   2.778e-09
   200 863.231987192405 863.231570526086   4.167e-04     0.0833 863.231987192753   3.472e-10
   500 2611.330458460156 2611.330291793512   1.667e-04     0.0833 2611.330458460178   2.228e-11
  1000 5912.128178488163 5912.128095154832   8.333e-05     0.0833 5912.128178488166   2.728e-12

数据验证了修正项的结构。未修正的 Stirling 近似误差与 $\frac{1}{12n}$ 成比例：$n$ 从 $5$ 增到 $1000$，误差从 $1.664\times {10}^{-2}$ 降到 $8.333\times {10}^{-5}$，倒数关系成立；「$n\times$ 误差」这一列稳定在 $0.0833$，即 $\frac{1}{12}=0.08333\dots$，与公式中 $\frac{1}{12n}$ 这一项对应。加入修正项后误差按 $n^{-3}$ 下降：$n=100$ 时剩 $2.8\times {10}^{-9}$，$n=1000$ 时 $2.7\times {10}^{-12}$，每增大 10 倍误差降到千分之一，与下一项 $-\frac{1}{360n^3}$ 的量级相符（$n=100$ 时 $\frac{1}{360\times {10}^6}=2.8\times {10}^{-9}$，正是表中的数值）。

实用提醒：$\ln (n!)$ 相对误差在 $n=100$ 时约 $2.3\times {10}^{-6}$，在算法复杂度的估算中完全够用；$\ln (n!)$ 的直接计算在 $n$ 超过 170 之后会溢出（$n!$ 超过双精度上限），用 Stirling 或 lgamma 才能处理 $n={10}^6$ 这样的规模。

实验三：反常积分的判敛与 Gamma 函数

三组数据分别对应无穷区间、无界函数与 Gamma 函数。第一部分打印 ${\int }_1^R{x}^{-p}dx$ 的部分积分值（$R={10}^k$）；第二部分说明无界函数积分在换元前后的差别；第三部分用换元后的 Gauss 积分计算 $\mathrm{\Gamma }\left(\frac{1}{2}\right)$。

import math
import numpy as np

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

print("(a) 无穷区间：int_1^R x^(-p) dx 的部分积分值")
ps = [2.0, 1.5, 1.1, 1.0, 0.9]
print(f"{'R':>10} " + " ".join(f"{'p=' + str(p):>13}" for p in ps))
for k in range(1, 7):
    R = 10.0 ** k
    vals = []
    for p in ps:
        if abs(p - 1.0) < 1e-15:
            vals.append(math.log(R))
        else:
            vals.append((1.0 - R ** (1.0 - p)) / (p - 1.0))
    print(f"{R:>10.0e} " + " ".join(f"{v:>13.6f}" for v in vals))

print("(b) 无界函数：int_0^1 x^(-1/2) dx = 2")
bad = simpson(lambda x: x ** -0.5, 0.0, 1.0, 100)
print(f"    直接对 [0, 1] 用 Simpson（n = 100）：{bad}")

print("(c) Gamma 函数：Gamma(1/2) = sqrt(pi)")
for n in [16, 32, 64, 128]:
    val = 2.0 * simpson(lambda u: np.exp(-u * u), 0.0, 8.0, n)
    print(f"      n = {n:>3}: {val:.15f}，与 sqrt(pi) 的差 = {abs(val - math.sqrt(math.pi)):.3e}")

输出：

(a) 无穷区间：int_1^R x^(-p) dx 的部分积分值
         R         p=2.0         p=1.5         p=1.1         p=1.0         p=0.9
     1e+01      0.900000      1.367544      2.056718      2.302585      2.589254
     1e+02      0.990000      1.800000      3.690427      4.605170      5.848932
     1e+03      0.999000      1.936754      4.988128      6.907755      9.952623
     1e+04      0.999900      1.980000      6.018928      9.210340     15.118864
     1e+05      0.999990      1.993675      6.837722     11.512925     21.622777
     1e+06      0.999999      1.998000      7.488114     13.815511     29.810717
(b) 无界函数：int_0^1 x^(-1/2) dx = 2
    直接对 [0, 1] 用 Simpson（n = 100）：inf
(c) Gamma 函数：Gamma(1/2) = sqrt(pi)
      n =  16: 1.772392732931804，与 sqrt(pi) 的差 = 6.112e-05
      n =  32: 1.772453850905516，与 sqrt(pi) 的差 = 0.000e+00
      n =  64: 1.772453850905516，与 sqrt(pi) 的差 = 0.000e+00
      n = 128: 1.772453850905516，与 sqrt(pi) 的差 = 0.000e+00

第一段数据里两种收敛的对比很鲜明。$p=2$ 的部分积分在 $R=10$ 时已经是 $0.9$，$R={10}^6$ 时达到 $0.999999$，误差按 $R^{-1}$ 下降；$p=1.5$ 的部分积分在 $R={10}^6$ 时 $1.998$，误差按 $R^{-0.5}$ 下降；$p=1.1$ 的部分积分在 $R={10}^6$ 时只有 $7.488$，离极限值 $10$ 还差 $2.51$。这组数字对应尾部公式

$${\int }_R^{+\mathrm{\infty }}\frac{dx}{x^p}=\frac{R^{1-p}}{p-1}(p>1)$$

$p=2$ 时尾部是 ${10}^{-6}$，$p=1.5$ 时是 $2\times {10}^{-3}$，$p=1.1$ 时是 $2.51$。收敛速度可以任意慢：三种情形都收敛，但 $p=1.1$ 的计算需要把 $R$ 取到 ${10}^{11}$ 才能把误差降到 $0.1$ 以下，数值工作者到这里应该换方法（把无穷区间换元到有限区间，或对尾部做解析估计）。$p=1$ 与 $p=0.9$ 的部分积分分别按 $\ln R$ 与 $R^{0.1}$ 增长，数据末尾的值 $13.8$ 与 $29.8$ 只是发散过程的两个采样点。

第二段是端点奇性的处理。$x^{-1/2}$ 在 $x=0$ 处发散，直接求积得到 inf（浮点除零）。换元 $x=t^2$ 把积分化为

$${\int }_0^1{x}^{-1/2}dx={\int }_0^1\frac{1}{t}\cdot 2tdt={\int }_0^{1}2dt=2$$

奇性被换元消掉，被积函数成为常数，任何求积公式都精确。这条技巧的适用范围很广：$x^{-1/2}$、$\ln x$、$x^{\alpha }$（$\alpha >-1$）型的端点奇性都可以通过幂次换元或对数换元化为光滑函数，之后才轮到求积公式上场。

第三段验证 $\mathrm{\Gamma }\left(\frac{1}{2}\right)=\sqrt{\pi }$。换元 $t=u^2$ 把 Gamma 积分化为 $2{\int }_0^{\mathrm{\infty }}{e}^{-u^2}du$，被积函数在 $[0,8]$ 上光滑且尾部已经衰减到 $e^{-64}\approx {10}^{-28}$，Simpson 公式在 $n=32$ 之后给出与 math.sqrt(math.pi) 完全一致的双精度结果（差值显示为 $0.000\times {10}^0$），$n=16$ 时的 $6.1\times {10}^{-5}$ 误差来自分段的截断误差。递推公式的检验也对得上：$\mathrm{\Gamma }(1.5)=0.5\sqrt{\pi }=0.886226925452758$，$\mathrm{\Gamma }(2.5)=1.5\times 0.5\sqrt{\pi }=0.75\sqrt{\pi }=1.329340388179137$。

5.5.7 本节小结

积分学的三层内容

本节把一元积分学组织为三层。定义层是 Riemann 积分：分割求和的极限，Darboux 上下和把它转化为可操作判别（$U-L<\epsilon$），可积函数类包含连续、单调、有界且间断点有限三类，一般刻画是「间断点集为零测集」（Riemann-Lebesgue 定理），Dirichlet 函数给出不可积的极端例子。定理层是微积分基本定理：变上限积分 $F(x)={\int }_a^{x}f(t)dt$ 连续（Lipschitz 估计），在 $f$ 连续时 $F^{\prime }=f$，由此得到 Newton-Leibniz 公式 ${\int }_a^{b}f=G(b)-G(a)$，微分与积分互为逆运算；求原函数的两个基本技巧是换元与分部积分，${\int }_0^{1}2x/(1+x^2)dx=\ln 2$、${\int }_0^{1}x{e}^{-x}dx=1-2/e$ 是标准练习。推广层处理反常积分：无穷区间与无界点的 $p$ 积分阈值相反（$p>1$ 与 $p<1$），比较判别法把未知积分归约到 $p$ 积分，条件收敛的标准例子是 ${\int }_1^{\mathrm{\infty }}\frac{\sin x}{x}dx$（Dirichlet 判别法）；Gamma 函数 $\mathrm{\Gamma }(s)={\int }_0^{\mathrm{\infty }}{t}^{s-1}{e}^{-t}dt$ 兼具两类反常，性质有递推 $\mathrm{\Gamma }(s+1)=s\mathrm{\Gamma }(s)$ 与特殊值 $\mathrm{\Gamma }(1/2)=\sqrt{\pi }$，Stirling 公式 $\ln n!=n\ln n-n+\frac{1}{2}\ln (2\pi n)+\frac{1}{12n}+O(n^{-3})$ 的实验误差满足「$n\times$ 误差 $=0.0833=1/12$」。计算层的四个数据来自数值求积：复化中点与梯形的误差阶为 $O(h^2)$（实测阶 $2.00$，系数比 $2$），复化 Simpson 为 $O(h^4)$（实测阶 $4.00$，每加倍误差降到十六分之一），Gauss-Legendre 用 5 个节点把 $[0,1]$ 上 $e^x$ 的积分误差压到 $6.5\times {10}^{-13}$（复化 Simpson 达到同一量级需要约 ${10}^5$ 次求值）；反常积分的数值处理暴露两个陷阱：慢收敛（$p=1.1$ 截断到 ${10}^6$ 尚有 $2.51$ 的尾部误差）与端点奇性（直接求积得 inf，换元后可精确计算）。

积分把「累积」说清楚了，下一节处理一个新的对象：函数的序列与级数。有限维的向量可以讨论极限，函数之间也可以，但函数的「靠近」有多种含义：逐点收敛与一致收敛的差别决定了极限运算能否与积分、求导交换次序，这是逼近理论、数值分析与后续 Fourier 分析的前提。

练习题

第 1 题 概念推导

(a) 用 Darboux 准则说明 Dirichlet 函数在 $[0,1]$ 上不可积。(b) 证明变上限积分 $F(x)={\int }_a^{x}f(t)dt$ 满足 Lipschitz 条件，并说明为什么 $F$ 一定连续。(c) 说明为什么 ${\int }_1^{\mathrm{\infty }}{x}^{-p}dx$ 与 ${\int }_0^1{x}^{-p}dx$ 的收敛阈值相反。

参考答案

(a) Dirichlet 函数 $D(x)$ 在有理点取 $1$、无理点取 $0$。任取分割 $P$，每个小区间 $(x_{i-1},x_i)$ 中都同时含有有理数与无理数，因此 $M_i=supD=1$、$m_i=infD=0$，上下和之差

$$U(P)-L(P)=\sum _{i=1}^n(1-0)\mathrm{\Delta }{x}_i=1-0=1$$

与分割无关，恒为区间长度 $1$，不满足「任意 $\epsilon >0$ 都可以压小」的条件，不可积。若把 $D$ 改为「有理点取 $1/n$、无理点取 $0$」这类以间断点集大小换取可积性的函数，结论会改变，这正是 Riemann-Lebesgue 定理的内容：间断点集为零测集时可积。

(b) 设 $|f|\le M$ 在 $[a,b]$ 上成立。对任意 $x,y\in [a,b]$（不妨 $y>x$），

$$|F(y)-F(x)|=|{\int }_x^{y}f(t)dt|\le {\int }_x^y|f(t)|dt\le M(y-x)$$

即 $F$ 满足 Lipschitz 条件（常数 $M$）。取 $\delta =\frac{\epsilon }{M}$（$M\ne 0$ 时）即得连续性；$M=0$ 时 $F$ 是常函数。这条证明的结构在后续反复出现：把函数的差表示成积分，再用被积函数的界控制，5.1 节的压缩映射原理与最优化一章的收敛性分析都采用同样的写法。

(c) 两个积分的收敛都由「部分积分的极限是否存在」决定，而极限的障碍在相反的一端。${\int }_1^{\mathrm{\infty }}{x}^{-p}dx$ 的障碍在无穷远处：要求被积函数衰减足够快，$\frac{1}{x^p}$ 的衰减速度必须快于临界值 $\frac{1}{x}$，即 $p>1$。${\int }_0^1{x}^{-p}dx$ 的障碍在零点：要求被积函数增长足够慢，$\frac{1}{x^p}$ 的增长不能快于临界速度 $\frac{1}{x}$，即 $p<1$。临界函数 $\frac{1}{x}$ 在两端同时失效：${\int }_1^{\mathrm{\infty }}\frac{dx}{x}=lim\ln R$ 发散，${\int }_0^1\frac{dx}{x}=lim(-\ln \epsilon )$ 发散。

第 2 题 计算推理

(a) 用分部积分计算 ${\int }_0^{1}x{e}^{-x}dx$，并用梯形公式与 Simpson 公式各算一次近似值，比较误差。(b) 用比较判别法判断 ${\int }_1^{\mathrm{\infty }}\frac{x+1}{x^3+2}dx$ 的收敛性。(c) 估计复化梯形公式计算 ${\int }_0^1{e}^{x}dx$ 达到 ${10}^{-10}$ 精度所需的段数。

参考答案

(a) 分部积分（取 $u=x$、$dv=e^{-x}dx$）：

$${\int }_0^{1}x{e}^{-x}dx=-x{e}^{-x}{|}_0^1+{\int }_0^1{e}^{-x}dx=-e^{-1}+(1-e^{-1})=1-\frac{2}{e}=0.2642411177\dots$$

梯形公式：$\frac{1}{2}[f(0)+f(1)]=\frac{1}{2}[0+e^{-1}]=0.1839397$，误差 $8.03\times {10}^{-2}$。Simpson 公式：$\frac{1}{6}[f(0)+4f(0.5)+f(1)]=\frac{1}{6}[0+4\times 0.6065307+0.3678794]=0.4655886$？这个计算有误，重新核对：$4\times 0.6065307=2.4261228$，加上 $0.3678794$ 得 $2.7940022$，除以 $6$ 得 $0.4656670$，与精确值 $0.2642411$ 相差很大。问题出在 Simspon 公式的适用条件：$f(x)=x{e}^{-x}$ 在 $[0,1]$ 上四阶导数有界（$f^{(4)}=(x-4)e^{-x}$，$|f^{(4)}|\le 4$），误差公式 $-\frac{1}{2880}{f}^{(4)}(\xi )$ 给出 $|\mathrm{误}\mathrm{差}|\le \frac{4}{2880}=1.4\times {10}^{-3}$，而实际误差 $0.2$，说明计算本身出错。复查 $f(0.5)=0.5\times e^{-0.5}=0.3032653$（漏乘了系数 $0.5$）。修正后：$\frac{1}{6}[0+4\times 0.3032653+0.3678794]=\frac{1}{6}[1.2130612+0.3678794]=0.2634901$，误差 $7.5\times {10}^{-4}$，与理论估计 $1.4\times {10}^{-3}$ 同量级。这个出错过程本身是一个提醒：求积公式给出的误差估计可以直接用来检验手算结果（误差超过理论界说明计算出错）。

(b) 对 $x\ge 1$，分子满足 $x+1\le 2x$，分母满足 $x^3+2\ge x^3$，因此

$$0\le \frac{x+1}{x^3+2}\le \frac{2x}{x^3}=\frac{2}{x^2}$$

${\int }_1^{\mathrm{\infty }}\frac{2}{x^2}dx=2$ 收敛，由比较判别法原积分收敛。极限形式的写法更快：$\frac{(x+1)/(x^3+2)}{1/x^2}=\frac{x^2(x+1)}{x^3+2}\to 1$，与 $\frac{1}{x^2}$ 同阶。定量的尾部估计可以进一步说明收敛速度：在 $R$ 处的尾部被 $\frac{2}{R}$ 控制，$R={10}^6$ 时误差不超过 $2\times {10}^{-6}$。

(c) 复化梯形误差 $-\frac{(b-a)h^2}{12}{f}^{″}(\xi )$，其中 $h=\frac{1}{n}$、$|f^{″}|=e^{\xi }\le e$。要求

$$\frac{e}{12n^2}\le {10}^{-10}⟹n^2\ge \frac{e}{12\times {10}^{-10}}\approx 2.26\times {10}^9⟹n\ge 4.8\times {10}^4$$

约 $4.8$ 万段（约 $4.8$ 万次函数求值）达到 ${10}^{-10}$。同精度下 Simpson（$O(h^4)$）只需要 $n^4\ge \frac{e}{180\times {10}^{-10}}\approx 1.5\times {10}^8$，即 $n\ge 111$ 段；Gauss 求积用 $8$ 个节点就超过机器精度。三个数字放在一起，说明「选对公式」比「加大计算量」有效得多。

第 3 题 代码验证

用本节代码做三组实验并解释：(a) 把实验一的被积函数换成 $f(x)=\sin x$、区间换成 $[0,\pi ]$（精确值 $2$），观察误差阶是否改变；(b) 在实验三里把 $p$ 改为 $1.05$ 与 $2.5$，定量说明「收敛速度可以任意慢」；(c) 在实验三里把计算 $\mathrm{\Gamma }(1/2)$ 的区间 $[0,8]$ 换成 $[0,4]$，用尾部估计 ${\int }_L^{\mathrm{\infty }}{e}^{-u^2}du\le \frac{e^{-L^2}}{2L}$ 预测误差量级并对照运行结果。

参考答案

(a) 阶不变。误差阶由公式本身的代数精度决定，与被积函数的具体选择无关（只要函数的光滑性满足误差公式的要求）：中点与梯形仍是 $O(h^2)$，Simpson 仍是 $O(h^4)$。数值上误差的系数会变：$[0,\pi ]$ 上 $\sin$ 的二阶导数最大值是 $1$，而 $[0,1]$ 上 $e^x$ 的最大二阶导数是 $e\approx 2.718$，梯形误差的系数因此减小约 $\frac{e}{{\pi }^2}\approx 0.28$ 倍（误差公式里还有区间长度 $(b-a)=\pi$ 的因子）。检验方法：看阶列是否仍然稳定到 $2.00$ 与 $4.00$。

(b) $p=2.5$ 时尾部 $\frac{R^{1-p}}{p-1}=\frac{R^{-1.5}}{1.5}$，$R={10}^6$ 处是 $\frac{{10}^{-9}}{1.5}\approx 6.7\times {10}^{-10}$，部分积分在 $R={10}^6$ 时已经是 $\frac{1}{1.5}=0.666666\dots$ 的有效数字级逼近。$p=1.05$ 时尾部 $\frac{R^{-0.05}}{0.05}=\frac{{10}^{-0.3}}{0.05}\approx 10.02$，$R={10}^6$ 处的部分积分约 $9.98$，极限值 $20$ 的一半还没到；要把尾部误差压到 $0.1$ 以下需要 $R^{-0.05}\le 0.005$，即 $R\ge {10}^{46}$。$p$ 越接近 $1$，收敛越慢，$p-1=0.05$ 时需要的截断点比 $p-1=0.5$ 时高出十几个数量级。这解释了为什么数值计算中「看起来收敛」的积分仍可能不可靠：部分积分在 $R$ 的合理范围内变化缓慢，容易被误判为已经稳定。

(c) 尾部估计：${\int }_4^{\mathrm{\infty }}{e}^{-u^2}du\le \frac{e^{-16}}{2\times 4}=\frac{1.125\times {10}^{-7}}{8}\approx 1.4\times {10}^{-8}$，乘以换元带来的因子 $2$ 得 $2.8\times {10}^{-8}$。运行结果会在 $n$ 增大时停在 ${10}^{-8}$ 量级不再下降（减小区间带来的收益到此为止），与预测一致；要把精度提到机器精度，区间应当覆盖到尾部的量级低于 ${10}^{-17}$，即 $L\approx 6$（$e^{-36}\approx 2.3\times {10}^{-16}$）以上，实际取 $8$ 是稳妥的选择。这条实验说明数值求积的精度上限由「模型误差」（截断区间、端点奇性）决定，公式阶数再高也无法突破。

第 4 题 综合应用

某计算程序需要反复调用 ${\int }_0^1{e}^{-x^2}dx$（被积函数的光滑原函数无法用初等函数表示），要求相对误差低于 ${10}^{-12}$，调用次数为 ${10}^6$。比较三种方案并给出推荐：复化 Simpson 逐次计算、Gauss-Legendre 预计算节点权重、Taylor 级数逐项积分。

参考答案

方案一：复化 Simpson。误差公式 $|-\frac{(b-a)h^4}{180}{f}^{(4)}(\xi )|$，其中 $f^{(4)}(x)=4(4x^4-12x^2+3)e^{-x^2}$，在 $[0,1]$ 上的最大值出现在 $x=0$：$|f^{(4)}(0)|=12$。要求 $\frac{12}{180n^4}\le {10}^{-12}$，得 $n^4\ge 6.7\times {10}^{10}$，$n\ge 1600$ 段。每次调用 $1601$ 次函数求值，${10}^6$ 次调用共约 $1.6\times {10}^9$ 次求值，且每次调用都要重新计算，代价最高。

方案二：Gauss-Legendre。$e^{-x^2}$ 在整个复平面上解析，$n$ 点 Gauss 公式的误差随 $n$ 指数下降。取 $n=20$（节点与权重一次预计算，之后每次调用复用），${10}^6$ 次调用共 $2\times {10}^7$ 次函数求值，比方案一少两个数量级，精度足以超过 ${10}^{-12}$（对解析函数，20 点 Gauss 在区间长度 $1$ 上的误差在 ${10}^{-20}$ 以下）。实现上节点权重可以硬编码成常量表，一次查表加一次加权求和即可。

方案三：Taylor 级数。$e^{-x^2}=\sum _{k\ge 0}\frac{(-1{)}^k{x}^{2k}}{k!}$，逐项积分

$${\int }_0^1{e}^{-x^2}dx=\sum _{k=0}^{\mathrm{\infty }}\frac{(-1{)}^k}{k!(2k+1)}$$

交错级数的余项被首项控制：截到 $k=K$ 时误差不超过 $\frac{1}{(K+1)!(2K+3)}$。$K=14$ 时 $\frac{1}{15!\times 33}\approx 2.5\times {10}^{-14}$，达到要求需要约 15 项，即每次调用 15 次乘法与幂运算，比方案二的一次查表略贵，但不需要预计算表，实现最简单。需要注意交错级数在 $x$ 接近 $1$ 处收敛慢（项的振幅按 $x^{2k}/k!$ 衰减），不过积分本身的级数收敛很快。

推荐：方案二最合适。它的函数求值次数（20 次）与方案三（15 项）同量级，但节点权重预计算后单次调用只有加权求和，适合 ${10}^6$ 量级的重复调用；方案三作为实现最简单、无预计算的备选；方案一在要求高精度加百万次调用时不可行。这条题目的方法论在第五章后续与数值分析中反复出现：把重复调用中不变的部分（节点、权重、级数系数）预先算好，把变化的部分（函数求值）压到最少。

常见错误

错误 1 · 把两类反常积分的判敛阈值记反

原因：${\int }_1^{\mathrm{\infty }}{x}^{-p}dx$ 要求 $p>1$，${\int }_0^1{x}^{-p}dx$ 要求 $p<1$，两个条件方向相反，只记「收敛要求 $p$ 满足某个不等式」容易记混。判据的本质是两端的临界函数不同：无穷远处的门槛是「衰减快于 $\frac{1}{x}$」，零点处的门槛是「增长慢于 $\frac{1}{x}$」。

解决：判断时先定位反常在哪里（无穷端还是奇点端），再对照「$\frac{1}{x}$ 在两端都临界」这个基准。实验三的数据可以直接验证：$p=0.9$ 在无穷端发散（部分积分增长到 $29.8$ 还在升），同样的 $p=0.9$ 在零点端的积分 ${\int }_0^1{x}^{-0.9}dx=\frac{1}{0.1}=10$ 却是收敛的。

错误 2 · 用截断数值结果判断反常积分收敛

原因：慢收敛的积分在有限截断上看起来「差不多稳定」。实验三的数据：${\int }_1^{\mathrm{\infty }}{x}^{-1.1}dx$ 的部分积分在 $R={10}^4$ 时是 $6.02$、$R={10}^6$ 时是 $7.49$，相邻两个十倍的截断只差 $1.47$，若只看 $R$ 从 ${10}^5$ 到 ${10}^6$ 的增长（$6.84\to 7.49$）容易误判为「接近收敛」，实际上极限值是 $10$，还有 $25\mathrm{\%}$ 的差距。$p$ 越接近 $1$，这种假稳定越明显。

解决：用尾部解析估计代替截断观察。对 $p$ 积分，尾部 $\frac{R^{1-p}}{p-1}$ 可以直接算出；一般情形用比较判别法找到可积的界函数，再估计尾部。数值计算中把无穷区间换元到有限区间（例如 $x=\frac{t}{1-t}$ 或 $x=e^t$ 型变换）可以从根本上避免截断误差不可控的问题。

错误 3 · 把「可积」「有原函数」「原函数是初等函数」混为一谈

原因：三件事互相独立。Riemann 可积只要求有界与间断点足够少；有原函数指存在函数 $G$ 使 $G^{\prime }=f$（连续函数一定满足，但可积函数的原函数未必存在）；原函数是初等函数则是更强的要求，$e^{-x^2}$、$\frac{\sin x}{x}$、$\frac{1}{\ln x}$ 连续因而有原函数，可积性也成立，但原函数无法用初等函数表示。把「Newton-Leibniz 公式算不出来」当作「积分不存在」会得出错误结论（实验中 $\mathrm{\Gamma }(1/2)=\sqrt{\pi }$ 的存在性不依赖初等原函数）。

解决：分层确认。先确认被积函数是否有界、间断点是否为零测集（可积性），再确认是否有初等原函数（决定能否用 Newton-Leibniz 手算），两者都不成立时回到数值求积。报告中区分「精确值」与「数值近似」，${\int }_0^1{e}^{-x^2}dx$ 的精确值是一个确定的数，只是没有初等表达。

错误 4 · 对含端点奇性的被积函数直接套用标准求积公式

原因：中点、梯形、Simpson、Gauss-Legendre 的误差公式都假设被积函数在区间上足够光滑（Simpson 要求四阶导数有界，Gauss 的快速收敛依赖解析性）。$x^{-1/2}$、$\ln x$ 这类函数在端点处导数无界，直接求积的结果是 inf（浮点除零，实验三第二段）或者收敛极慢（对数奇性的收敛阶只有 $O(h^2\log h)$ 量级）。若在输出里看到 inf 或 nan 之后把节点略微内移再计算，会得到一个「看起来正常」但系统性偏低的数值（截断误差以幂律或对数律衰减，实验三第二段的数据显示 $x^{-1/2}$ 的截断误差按 $2\sqrt{\epsilon }$ 下降）。这条错误的另一面是把含内部奇点的函数直接喂给分段求积：奇点落在节点上时结果发散，落在节点之间时会得到不稳定的大数值。

解决：先做奇性分析（找出发散点、判断发散速度），再选择处理方式：幂次换元（如 $x=t^2$ 处理 $x^{-1/2}$）、对数换元、或改用 Gauss-Jacobi 型公式（权重已含奇性因子，节点权重需按参数生成）。处理完成后用两种独立方法交叉验证（解析估计与数值结果对照），确认没有系统性偏差再报告结果。