5.12 渐近分析与特殊函数

前面的内容处理的是精确问题：求极限、求积分、求极值。实际工作中另一类问题出现得更频繁：参数很大或很小时，精确解无法求得，但解的主项可以写成简单的显式表达式。例如 ${\int }_0^1{e}^{-n{x}^2}dx$ 对任意 $n$ 都有精确值 $\frac{\sqrt{\pi }}{2\sqrt{n}}\mathrm{erf}(\sqrt{n})$，但当 $n={10}^6$ 时，更有用的是 $\frac{\sqrt{\pi }}{2\sqrt{n}}$ 这个近似（相对误差 ${10}^{-400000}$ 量级的修正），它一眼就能看出「积分随 $n^{-1/2}$ 衰减」。渐近分析的目标就是得到这种表达式：不要精确解，要主项、修正项与误差量级。

三条线索构成本节。Laplace 方法处理形如 $\int e^{nf(x)}dx$ 的积分，$n$ 很大时积分值几乎完全由 $f$ 的最大值附近的区域决定，主项是 $e^{nf(x_0)}\sqrt{2\pi /(n|f^{″}(x_0)|)}$，修正项按 $1/n$ 的幂展开。Euler-Maclaurin 公式把求和化成积分加修正项，是分析级数与求和量的标准工具，调和数 $H_N=\ln N+\gamma +\frac{1}{2N}-\frac{1}{12N^2}+\cdots$ 是它的典型应用。奇异摄动与边界层处理「最高阶导数乘以小参数」的微分方程：解在大部分区域由低阶方程描述（外解），但在某个边界附近的宽度 $O(\epsilon )$ 的薄层内发生急剧变化（内解），两者用匹配条件拼成一致有效的复合展开。实验四给出的数据（外解误差 $1.5$ 量级、复合展开误差降到 ${10}^{-16}$）说明了这条方法的威力，同时揭示了它的数值困难：当网格步长大于 $2\epsilon$ 时，常规的中心差分格式会产生非物理振荡（实验四实测 $6$ 次斜率变号），必须改用指数拟合等保持单调的格式。

5.12.1 Laplace 方法

主项的来源

考虑积分

$$I(n)={\int }_a^b{e}^{nf(x)}dx$$

$n\to +\mathrm{\infty }$ 时，被积函数在 $f$ 取最大值的点 $x_0$ 附近有一个宽度 $O(n^{-1/2})$ 的尖峰，其余部分的贡献相对可以忽略。把 $f$ 在 $x_0$ 处展开（设 $x_0$ 是内部点、$f^{″}(x_0)<0$）

$$f(x)=f(x_0)+\frac{1}{2}{f}^{″}(x_0)(x-x_0{)}^2+O((x-x_0{)}^3)$$

则

$$I(n)\approx e^{nf(x_0)}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{\frac{n}{2}{f}^{″}(x_0)(x-x_0{)}^2}dx=e^{nf(x_0)}\sqrt{\frac{2\pi }{n|f^{″}(x_0)|}}$$

把积分区间从 $[a,b]$ 换成 $\mathbb{R}$ 的代价是 $O(e^{-cn})$ 量级的（高斯尾部），这条误差比任何 $n$ 的幂次都小，因此在幂次展开的语境下可以忽略。三条常见变体：

端点最大值：$f$ 在端点 $a$ 取最大且 $f^{\prime }(a)<0$，则 $I(n)\approx \frac{e^{nf(a)}}{n|f^{\prime }(a)|}$，衰减变成 $1/n$ 而不是 $n^{-1/2}$（换元 $t=n(x-a)$ 后积分到无穷，被积函数指数衰减，量级是 $1/n$）。

多个最大值：$f$ 在若干点取到同一个最大值时，主项是各点贡献之和；这些点的贡献量级相同，不能只保留一个。

最大值处 $f^{″}(x_0)=0$：此时的峰更平坦，需要展到更高阶。设 $f(x)-f(x_0)\sim -c|x-x_0{|}^{2m}$，主项是 $O(n^{-1/(2m)})$（用 $\mathrm{\Gamma }$ 函数表达）。实验一第一组用的 ${\int }_0^1{e}^{-n{x}^2}dx$ 属于端点最大值情形（最大值在 $x=0$ 且 $f^{\prime }=0$，属于半高斯），主项 $\frac{\sqrt{\pi }}{2\sqrt{n}}$；第二组的 ${\int }_{-\pi }^{\pi }{e}^{a\cos x}dx$ 属于内部最大值情形。

修正项

把 $f$ 展开到四阶、并把 $e^{n\cdot O((x-x_0{)}^4)}$ 幂级数展开，可以算出修正项。对 $f(x)=\cos x$（$x_0=0$、$f^{″}(0)=-1$、$f^{(4)}(0)=1$），得到

$${\int }_{-\pi }^{\pi }{e}^{a\cos x}dx=e^a\sqrt{\frac{2\pi }{a}}(1+\frac{1}{8a}+\frac{9}{128a^2}+\cdots )$$

实验一第二组的数据验证这条展开：$a=1$ 时主项的相对误差 $1.435\times {10}^{-1}$，加 $1/(8a)$ 后降到 $3.639\times {10}^{-2}$；$a=10$ 时从 $1.312\times {10}^{-2}$ 降到 $7.803\times {10}^{-4}$；$a=20$ 时从 $6.395\times {10}^{-3}$ 降到 $1.845\times {10}^{-4}$。主项的误差按 $1/a$ 下降，加一项修正后按 $1/a^2$ 下降（$a$ 从 $10$ 增到 $20$，$7.803\times {10}^{-4}$ 与 $1.845\times {10}^{-4}$ 的比值是 $4.23$，接近 $2^2=4$）。这条规律是渐近展开的普遍性质：每一项把误差的幂次提高一次（或两次，取决于展开的对称性）。

第一组的数据给出指数级精度的例子：$I(n)={\int }_0^1{e}^{-n{x}^2}dx$ 的主项误差乘上 $e^n\sqrt{\pi n}$ 后趋于 $1$（$n=5$ 时 $0.9222$，$n=10$ 时 $0.9561$，$n=20$ 时 $0.9767$），说明误差的主项确实是 $\frac{e^{-n}}{\sqrt{\pi n}}$。$n=20$ 时主项的相对误差是 $2.54\times {10}^{-10}$，已经低于一般数值积分的精度——这是渐近展开相对于数值积分的优势：参数大时渐近公式比数值方法更精确，代价是它只对特定的参数范围有效。

鞍点法与最陡下降

把积分推广到复平面上的围道积分 ${\int }_C{e}^{n\varphi (z)}dz$，最大值原理变成鞍点原理：被积函数的模在鞍点处沿最陡下降方向取极大，把围道变形到经过鞍点的最陡下降路径上，再作 Laplace 方法（复版本）。结论是

$${\int }_C{e}^{n\varphi (z)}dz\sim \frac{e^{n\varphi (z_0)}\sqrt{2\pi }}{\sqrt{-n{\varphi }^{″}(z_0)}}\cdot (\text{符号因子})$$

鞍点法的应用范围超出纯数学：统计力学中的配分函数、组合计数中的系数提取（生成函数的系数 $[x^n]F(x)=\frac{1}{2\pi i}\oint F(z)z^{-n-1}dz$，$n$ 大时用鞍点法）、概率论中的大偏差原理（第四章 4.6 节的速率函数由 Legendre 变换给出，与鞍点方程 ${\varphi }^{\prime }(z_0)=0$ 是同一个条件）都是这条方法的实例。$\mathrm{\Gamma }$ 函数的 Stirling 公式、中心二项式系数 $(\genfrac{}{}{0ex}{}{2n}{n})\sim \frac{4^n}{\sqrt{\pi n}}$、熵与组合数的指数率都可以由鞍点法统一导出。

5.12.2 Euler-Maclaurin 公式

公式

Euler-Maclaurin 公式把求和与积分联系起来：

$$\sum _{k=1}^{N}f(k)={\int }_1^{N}f(x)dx+\frac{f(1)+f(N)}{2}+\sum _{j=1}^m\frac{B_{2j}}{(2j)!}(f^{(2j-1)}(N)-f^{(2j-1)}(1))+R_m$$

其中 $B_{2j}$ 是 Bernoulli 数（$B_2=\frac{1}{6}$、$B_4=-\frac{1}{30}$、$B_6=\frac{1}{42}$、$\dots$），$R_m$ 是余项。公式的含义是：求和可以用积分加上端点修正与导数修正来逼近，修正项的系数由 Bernoulli 数给出。对 $f(x)=\frac{1}{x}$ 用这条公式（把下限换成 $1$、注意 $f$ 的高阶导数），得到调和数的展开

$$H_N=\sum _{k=1}^N\frac{1}{k}=\ln N+\gamma +\frac{1}{2N}-\frac{1}{12N^2}+\frac{1}{120N^4}-\frac{1}{252N^6}+\cdots$$

实验二的数据给出四项的完整误差阶梯：只用 $\ln N+\gamma$ 时误差是 $4.917\times {10}^{-2}$（$N=10$）、$4.992\times {10}^{-3}$（$N=100$）、$5.000\times {10}^{-7}$（$N={10}^6$），严格按 $O(1/N)$ 下降（系数 $\frac{1}{2}$）；加上 $\frac{1}{2N}$ 后误差变成 $8.325\times {10}^{-4}$、$8.333\times {10}^{-6}$、$8.349\times {10}^{-14}$，按 $O(1/N^2)$ 下降且系数等于 $\frac{1}{12}$（$N=100$ 时 $\frac{1}{12\times {10}^4}=8.33\times {10}^{-6}$ ✓）；再减 $\frac{1}{12N^2}$ 后误差按 $O(1/N^4)$ 下降，系数 $\frac{1}{120}$（$N=10$ 时 $8.294\times {10}^{-7}$ 对 $\frac{1}{120\times {10}^4}=8.33\times {10}^{-7}$ ✓）；加 $\frac{1}{120N^4}$ 后按 $O(1/N^6)$ 下降，系数 $\frac{1}{252}$（$N=10$ 时 $3.927\times {10}^{-9}$ 对 $\frac{1}{252\times {10}^6}=3.97\times {10}^{-9}$ ✓）。

这组数据把「渐近展开」的准确含义写清楚了：每一项的系数可以精确预测，误差与下一项同阶。同时要注意展开式的发散性：Bernoulli 数 $B_{2j}$ 随 $j$ 增长极快（$|B_{2j}|\sim 2(2j)!/(2\pi {)}^{2j}$），级数对固定的 $N$ 最终会发散；实用的做法是取到「误差最小」的那一项为止（截断误差与舍入误差的平衡，与 5.4 节数值微分的最优步长是同一类权衡）。这一点与 Taylor 级数的绝对收敛形成对照：Taylor 级数收敛但可能收敛得很慢，渐近级数发散但在参数大时极其精确。

应用

三条常见用法。

级数求和：把 $\sum f(k)$ 化成积分加修正，用于数值求和的加速（实验二的做法），也用于估计级数尾部的量级（5.2 节的比较判别法给出方向，Euler-Maclaurin 给出常数）。

渐近展开的推导：Stirling 公式的对数形式 $\ln n!=n\ln n-n+\frac{1}{2}\ln (2\pi n)+\frac{1}{12n}-\frac{1}{360n^3}+\cdots$ 就是 Euler-Maclaurin 公式作用于 $\ln n!=\sum _{k=1}^n\ln k$ 的结果（5.5 节给出了前两项的数值验证）。

离散与连续的桥：Euler-Maclaurin 公式是「离散和 ↔ 连续积分」的定量版本，在数值积分（梯形公式的误差展开）、数论（$\sum _{k\le n}d(k)$ 的估计）、算法分析（调和级数型的代价估计）中反复出现。

5.12.3 特殊函数的渐近与递推

三类特殊函数

Gamma 函数 $\mathrm{\Gamma }(s)={\int }_0^{\mathrm{\infty }}{t}^{s-1}{e}^{-t}dt$（5.5 节）：用 Laplace 方法（被积函数的对数在 $t=s-1$ 处取极大）得到 Stirling 公式，相对误差按 $1/s$ 的幂下降。

误差函数 $\mathrm{erf}(z)=\frac{2}{\sqrt{\pi }}{\int }_0^z{e}^{-t^2}dt$ 与余误差函数 $\mathrm{erfc}(z)=1-\mathrm{erf}(z)$：$z$ 大时 $\mathrm{erfc}(z)$ 极小，用 Laplace 方法（或分部积分）得

$$\mathrm{erfc}(z)\sim \frac{e^{-z^2}}{z\sqrt{\pi }}(1-\frac{1}{2z^2}+\frac{3}{4z^4}-\cdots )$$

这条展开在统计中直接对应「正态分布尾概率」的估计（第二章的尾概率与第六章的假设检验），也是大偏差原理的具体形态：$z$ 大时尾概率按 $e^{-z^2}$ 指数衰减，前面的 $1/z$ 因子与级数修正是幂次修正。

Bessel 函数 $J_{\nu }(x)$：柱坐标下波动方程与热传导方程的解，$x$ 大时

$$J_{\nu }(x)\sim \sqrt{\frac{2}{\pi x}}\cos (x-\frac{\nu \pi }{2}-\frac{\pi }{4})$$

即在远处表现为振幅按 $x^{-1/2}$ 衰减的振荡。$x$ 小时 $J_{\nu }(x)\sim \frac{1}{\mathrm{\Gamma }(\nu +1)}{\left(\frac{x}{2}\right)}^{\nu }$。

递推关系的稳定性

Bessel 函数满足三项递推

$$J_{n+1}(x)=\frac{2n}{x}{J}_n(x)-J_{n-1}(x)$$

它在理论上对一切整数 $n$ 成立，但数值上不能随意使用。原因在于递推有两个解：$J_n$（随 $n$ 增大迅速衰减）与 $Y_n$（随 $n$ 增大迅速增长）。用前向递推从 $J_0$、$J_1$ 出发时，任何舍入误差都会被 $Y_n$ 的分量放大（放大因子约 ${\left(\frac{2n}{x}\right)}^n$），几步之后结果完全被误差淹没。

实验三的数据（$x=1$、$n$ 到 $20$）给出完整的量化形态：前向递推的相对误差从 $n=5$ 时的 $8.03\times {10}^{-12}$（还算可用）增到 $n=8$ 时的 $3.48\times {10}^{-5}$、$n=10$ 时的 $3.561$（已经完全没有意义）、$n=12$ 时的 $8.20\times {10}^5$、$n=20$ 时的 $8.18\times {10}^{29}$。后向递推（Miller 算法）从 $n=60$ 处的一个任意小值出发向下递推，再用恒等式 $J_0(x)+2\sum _{k\ge 1}{J}_{2k}(x)=1$ 归一化，得到的相对误差始终在 ${10}^{-15}$ 量级（$n=20$ 时 $1.897\times {10}^{-15}$）。

这条对照给出常微分方程与递推计算的一般原则：递推的方向应当与解的衰减方向一致（向前递推用于增长解，向后递推用于衰减解）。同一条原则出现在很多场合：三对角矩阵的追赶法（向前消元对应稳定的方向）、卡尔曼滤波的协方差更新（用 Joseph 形式保持正定）、偏微分方程的迎风格式（实验四）。判断稳定性的方法是检查递推的放大因子与解的增长/衰减趋势是否匹配，而不是只看公式是否成立。

5.12.4 奇异摄动与边界层

问题的结构

考虑含小参数的微分方程

$$\epsilon y^{″}+y^{\prime }=f(x),y(0)=\alpha ,y(1)=\beta$$

$\epsilon \to 0$ 时方程「降阶」为一阶方程，而两个边界条件仍然存在，这个矛盾说明解在某个边界附近必然出现急剧变化。以实验四的例子（$f(x)=x$、$\alpha =0$、$\beta =2$）说明。

外解：在远离边界的区域，忽略 $\epsilon y^{″}$ 项，解 $y_1^{\prime }=x$ 得 $y_1=\frac{x^2}{2}+C$。用哪个边界条件定 $C$？取远处的边界 $x=1$（这里假设层在 $x=0$）：$y_1=\frac{x^2}{2}+\frac{3}{2}$。

内解：在 $x=0$ 附近用放大坐标 $\xi =\frac{x}{\epsilon }$，方程变成 $\frac{d^{2}Y}{d{\xi }^2}+\frac{dY}{d\xi }=0$（$f$ 的项被放大后量级更高），解为 $Y=A+B{e}^{-\xi }$。边界条件 $y(0)=0$ 给出 $A+B=0$。

匹配：内解在 $\xi \to \mathrm{\infty }$ 时的极限要与外解在 $x\to 0$ 时的极限一致（匹配条件）：$A=y_1(0)=\frac{3}{2}$，因此 $B=-\frac{3}{2}$。

复合展开：把外解与内解相加再减去公共部分（匹配值），得一致有效的近似

$$y_{\text{comp}}(x)=\underset{y_{\text{out}}}{\underbrace{\frac{x^2}{2}+\frac{3}{2}+\epsilon (1-x)}}-\frac{3}{2}{e}^{-x/\epsilon }$$

（外解中保留一阶修正 $\epsilon (1-x)$ 是标准做法，它使误差下降到指数级。）

实验四的数据验证这条构造：只用外解时最大误差是 $1.600$（$\epsilon =0.1$）、$1.510$（$\epsilon =0.01$），把「层宽趋于零、误差不趋于零」这一点显示得很清楚——误差集中在 $x=0$ 附近，正是被忽略的边界层。复合展开的最大误差是 $7.264\times {10}^{-5}$（$\epsilon =0.1$）、$3.195\times {10}^{-9}$（$\epsilon =0.05$）、$4.441\times {10}^{-16}$（$\epsilon \le 0.02$，机器精度），误差随 $\epsilon$ 减小按指数下降（$\epsilon$ 减半时误差降到约 $\frac{1}{{10}^4}$），与「剩余的误差只有 $e^{-\text{常数}/\epsilon }$ 量级」这条理论预期一致。

边界层与数值振荡

渐近分析解决了「解的形状」，但数值求解仍然困难：边界层的宽度 $O(\epsilon )$ 意味着在均匀网格上要么把步长取到 $O(\epsilon )$（代价可能是 ${10}^4$ 个网格点以上），要么承受解的急剧变化带来的格式问题。

实验四的数值部分给出三种格式的完整对比（$\epsilon =0.01$）：

中心差分：方程为 $\frac{\epsilon (y_{i+1}-2y_i+y_{i-1})}{h^2}+\frac{y_{i+1}-y_{i-1}}{2h}=x_i$。网格 Péclet 数 $Pe=\frac{h}{2\epsilon }$ 大于 $1$（即 $h>2\epsilon$）时，一阶项的系数 $-\frac{1}{2h}$ 使三对角矩阵出现负的非对角元，格式失去单调性，解出现振荡。$N=20$（$h=0.05$、$h/(2\epsilon )=2.5$）时最大误差 $0.657317$，解中有 $6$ 处斜率变号（非物理振荡）；$N=50$ 及更细的网格上（$h/(2\epsilon )\le 1$）斜率变号消失，误差降到 $0.204356$ 及以下。

迎风格式：一阶项改用后向差分（对流方向来自左侧），系数变为 $\frac{\epsilon }{h^2}-\frac{1}{h}$，在 $h<2\epsilon$ 时为正、$h=2\epsilon$ 时为零、$h>2\epsilon$ 时为负。当 $h=2\epsilon$（$N=50$）时对角元为零，矩阵奇异，方程无解——这是实验中唯一一次求解失败，也是稳定性条件在数值上的直接体现。$h<2\epsilon$ 时格式单调但只有一阶精度（$N=2000$ 时误差 $0.014026$，比中心差分的 $0.000116$ 差两个数量级）。

指数拟合格式（Scharfetter-Gummel）：利用单元内方程 $\epsilon y^{″}+y^{\prime }=f$ 的精确解（对 $f$ 取单元内的常值近似）导出系数

$$\frac{\epsilon }{h^2}[y_{i+1}-(1+\beta )y_i+\beta y_{i-1}]=f_i\frac{\epsilon }{h}(1-\beta ),\beta =e^{-h/\epsilon }$$

三个系数 $\frac{\epsilon \beta }{h^2}$、$-\frac{\epsilon (1+\beta )}{h^2}$、$\frac{\epsilon }{h^2}$ 恒为正（或负对角），格式对任意 $h$ 保持单调，误差在全部网格上最小：$N=20$ 时 $0.014469$、$N=100$ 时 $0.000773$、$N=2000$ 时 $0.000002$，且没有任何振荡。代价是系数需要计算指数函数（对极小 $\beta =e^{-h/\epsilon }$ 要注意下溢，可以在 $\beta$ 小于某个阈值时改用零阶近似）。

这组数据的方法论含义超出这个方程：当解在某个区域急剧变化时，格式的单调性与稳定性比形式上的精度阶更重要。中心差分的二阶精度只在 $h<2\epsilon$ 时有效，而这个条件在 $\epsilon ={10}^{-4}$ 时要求 $N>50000$；指数拟合格式在粗网格上就能给出可用结果，代价是精度只有一阶到二阶（在 $h$ 与 $\epsilon$ 相当时）。工程实践中常把两种手段结合：用指数拟合或迎风保证全局稳定，在感兴趣的区域内加密网格或改用高阶格式。

5.12.5 渐近分析在算法与统计中的用法

三条常见的应用形态。

算法复杂度：Euler-Maclaurin 公式给出求和量的精确展开，例如二叉堆建堆的代价 $\sum _{h=0}^{\log n}\frac{n}{2^{h+1}}\cdot O(h)$ 的估计、快速排序的平均比较次数 $2n\ln n+O(n)$ 中常数项的确定、哈希表长链的概率估计。复杂度分析需要「主项 + 低阶项」这一对信息，渐近展开正好提供。

大偏差与尾概率：$\int e^{nf(x)}dx$ 的 Laplace 方法与概率中的大偏差原理是同一回事（第四章 4.6 节的 Cramér 定理）：把概率写成积分的指数形式，用 Laplace 方法提出主项，得到 $P({\overline{X}}_n\ge a)\approx e^{-nI(a)}$ 的形式，速率函数 $I(a)=\underset{t}{sup}\{ta-\ln M(t)\}$ 由 Legendre 变换给出（$M$ 是矩生成函数）。这条对应关系说明「积分的渐近」与「概率的指数衰减」在数学上是同一个机制。

统计中的近似：Stirling 公式用于二项系数的估计（第四章的组合数与熵），$\mathrm{erfc}$ 的展开用于正态尾概率，Beta 与 Gamma 函数的比值用于贝叶斯因子的近似（第六章）。这些近似的共同前提是「某个参数大」（样本量、维度、信噪比），使用时要说明参数的取值范围与误差量级，因为渐近公式在小参数处可能完全失效（实验一第二组在 $a=1$ 时主项的相对误差是 $14\mathrm{\%}$）。

5.12.6 数值实验

实验一：Laplace 方法的主项与修正项

两组积分分别对应端点最大值与内部最大值的情形，用精确值核对主项与第一修正项。

import math
from scipy.special import erf, i0

print("(a) I(n) = int_0^1 exp(-n x^2) dx，精确值 = sqrt(pi)/(2 sqrt(n)) * erf(sqrt(n))")
print(f"{'n':>6} {'精确值':>16} {'主项 sqrt(pi/n)/2':>20} {'相对误差':>14} "
      f"{'误差 * e^n * sqrt(pi n)':>26}")
for n in [1, 2, 3, 5, 10, 20]:
    exact = math.sqrt(math.pi) / (2 * math.sqrt(n)) * erf(math.sqrt(n))
    lead = math.sqrt(math.pi) / (2 * math.sqrt(n))
    rel = abs(lead - exact) / exact
    print(f"{n:>6} {exact:>16.10f} {lead:>20.10f} {rel:>14.3e} "
          f"{rel * math.exp(n) * math.sqrt(math.pi * n):>26.6f}")

print("(b) I(a) = int_{-pi}^{pi} exp(a cos x) dx = 2 pi I_0(a)，最大值在 x = 0")
print("    Laplace 展开：I(a) = e^a sqrt(2 pi / a) (1 + 1/(8a) + 9/(128 a^2) + ...)")
print(f"{'a':>6} {'精确值 2 pi I_0(a)':>22} {'主项':>20} {'主项相对误差':>16} {'加 1/(8a) 后':>16}")
for a in [1, 2, 5, 10, 20]:
    exact = 2 * math.pi * float(i0(a))
    lead = math.exp(a) * math.sqrt(2 * math.pi / a)
    lead2 = lead * (1 + 1 / (8 * a))
    print(f"{a:>6} {exact:>22.6f} {lead:>20.6f} {abs(lead - exact) / exact:>16.3e} "
          f"{abs(lead2 - exact) / exact:>16.3e}")

输出：

(a) I(n) = int_0^1 exp(-n x^2) dx，精确值 = sqrt(pi)/(2 sqrt(n)) * erf(sqrt(n))
     n              精确值      主项 sqrt(pi/n)/2           相对误差      误差 * e^n * sqrt(pi n)
     1     0.7468241328         0.8862269255      1.867e-01                   0.899337
     2     0.5981440067         0.6266570687      4.767e-02                   0.882911
     3     0.5043435602         0.5116633540      1.451e-02                   0.894935
     5     0.3957123096         0.3963327298      1.568e-03                   0.922229
    10     0.2802473905         0.2802495608      7.744e-06                   0.956094
    20     0.1981663648         0.1981663649      2.540e-10                   0.976674
(b) I(a) = int_{-pi}^{pi} exp(a cos x) dx = 2 pi I_0(a)，最大值在 x = 0
    Laplace 展开：I(a) = e^a sqrt(2 pi / a) (1 + 1/(8a) + 9/(128 a^2) + ...)
     a        精确值 2 pi I_0(a)                   主项           主项相对误差       加 1/(8a) 后
     1               7.954927             6.813722        1.435e-01        3.639e-02
     2              14.323057            13.096761        8.562e-02        2.847e-02
     5             171.153162           166.370891        2.794e-02        3.640e-03
    10           17691.669349         17459.618631        1.312e-02        7.803e-04
    20       273684760.984164     271934665.877386        6.395e-03        1.845e-04

(a) 端点最大值的指数级精度。$n=20$ 时主项已经与精确值在十位小数上一致（$0.1981663648$ 对 $0.1981663649$，相对误差 $2.54\times {10}^{-10}$），$n=1$ 时的相对误差是 $1.867\times {10}^{-1}$。最后一列（误差乘以 $e^n\sqrt{\pi n}$）的值从 $0.899337$ 缓慢升到 $0.976674$，说明误差的主项确实是 $\frac{e^{-n}}{\sqrt{\pi n}}$，剩余偏差来自 $1/n$ 的幂次修正。这条数据的实用含义是渐近公式的精度随参数以指数速度提升：$n$ 从 $1$ 增到 $20$，相对误差从 $18.7\mathrm{\%}$ 降到 $2.5\times {10}^{-10}$，跨越了九个数量级。

(b) 内部最大值的幂次修正。主项的相对误差按 $1/a$ 下降（$a=1$ 时 $1.435\times {10}^{-1}$、$a=20$ 时 $6.395\times {10}^{-3}$，比值约 $22$ 对应 $a$ 的比值 $20$），加 $1/(8a)$ 修正后按 $1/a^2$ 下降（$3.639\times {10}^{-2}$ 到 $1.845\times {10}^{-4}$，比值 $197$ 对应 $a$ 的平方比 $400$ 的量级）。小参数处的失效也很清楚：$a=1$ 时即使加了修正项，相对误差仍有 $3.6\mathrm{\%}$，因为 $a$ 不够大，展开的参数（$1/a$）不满足小量条件。两条数据合起来说明渐近展开的适用条件：参数必须落在展开参数为小量的范围内，超出范围时应当改用数值积分或保留更多项（如果收敛性允许）。

实验二：Euler-Maclaurin 公式求和

用四项修正计算调和数 $H_N$，观察误差幂次逐级提高。

import math

print("H_N = sum_{k=1}^N 1/k = ln N + gamma + 1/(2N) - 1/(12 N^2) + 1/(120 N^4) - ...")
gamma = 0.5772156649015329
rows = []
for N in [10, 100, 1000, 10 ** 4, 10 ** 6]:
    s = math.fsum(1.0 / k for k in range(1, N + 1))
    lN = math.log(N)
    v1 = lN + gamma
    v2 = v1 + 1.0 / (2 * N)
    v3 = v2 - 1.0 / (12 * N * N)
    v4 = v3 + 1.0 / (120 * N ** 4)
    rows.append((N, s, v1, v2, v3, v4))
print(f"{'N':>8} {'|ln N + gamma|':>18} {'加 1/(2N)':>16} {'减 1/(12N^2)':>18} "
      f"{'加 1/(120N^4)':>18}")
for N, s, v1, v2, v3, v4 in rows:
    print(f"{N:>8} {abs(v1 - s):>18.3e} {abs(v2 - s):>16.3e} {abs(v3 - s):>18.3e} "
          f"{abs(v4 - s):>18.3e}")

输出：

H_N = sum_{k=1}^N 1/k = ln N + gamma + 1/(2N) - 1/(12 N^2) + 1/(120 N^4) - ...
       N     |ln N + gamma|         加 1/(2N)        减 1/(12N^2)       加 1/(120N^4)
      10          4.917e-02        8.325e-04          8.294e-07          3.927e-09
     100          4.992e-03        8.333e-06          8.333e-11          3.553e-15
    1000          4.999e-04        8.333e-08          8.882e-15          8.882e-16
   10000          5.000e-05        8.333e-10          1.776e-15          1.776e-15
 1000000          5.000e-07        8.349e-14          0.000e+00          0.000e+00

四列的误差衰减速度依次是 $O(1/N)$、$O(1/N^2)$、$O(1/N^4)$、$O(1/N^6)$。第一列的数值可以直接核对系数：$N=10$ 时 $4.917\times {10}^{-2}$，而 $\frac{1}{2N}=0.05$，两者相符（误差的主项就是被丢掉的第一项）；$N={10}^6$ 时 $5.000\times {10}^{-7}$，正是 $\frac{1}{2\times {10}^6}$。第二列的数值等于 $\frac{1}{12N^2}$：$N=10$ 时 $8.325\times {10}^{-4}$ 对 $\frac{1}{1200}=8.333\times {10}^{-4}$，$N=100$ 时 $8.333\times {10}^{-6}$ 对 $\frac{1}{12\times {10}^4}$ 精确相符。第三列对应 $\frac{1}{120N^4}$（$N=10$ 时 $8.294\times {10}^{-7}$ 对 $\frac{1}{120\times {10}^4}=8.333\times {10}^{-7}$），第四列对应下一项 $\frac{1}{252N^6}$（$N=10$ 时 $3.927\times {10}^{-9}$ 对 $\frac{1}{252\times {10}^6}=3.97\times {10}^{-9}$）。

$N=1000$ 之后第三、四列落到 ${10}^{-15}$ 量级并保持不变，这是双精度舍入的下限（$H_N$ 本身用 math.fsum 精确求和，误差只剩舍入误差）。这条现象给出渐近展开的实用边界：展开的项数应当取到误差与舍入误差相当为止，继续增加项数没有收益（对发散的渐近级数还会使误差回升）。

实验三：特殊函数的递推稳定性

用前向与后向递推计算 Bessel 函数 $J_n(1)$，与库函数对照。

from scipy.special import jv

x = 1.0
NMAX = 20
j0 = float(jv(0, x))
j1 = float(jv(1, x))
fw = [j0, j1]
for n in range(1, NMAX):
    fw.append((2 * n / x) * fw[n] - fw[n - 1])       # 前向递推
exact = [float(jv(n, x)) for n in range(NMAX + 1)]
Nstart = 60
b = [0.0] * (Nstart + 2)
b[Nstart] = 1e-30
for n in range(Nstart, 0, -1):
    b[n - 1] = (2 * n / x) * b[n] - b[n + 1]         # 后向递推（Miller 算法）
norm = b[0] + 2 * sum(b[2::2])                       # 归一化：J_0 + 2 sum J_2k = 1
bw = [v / norm for v in b]
print(f"{'n':>4} {'精确值':>16} {'前向递推':>18} {'前向相对误差':>16} "
      f"{'后向递推':>16} {'后向相对误差':>16}")
for n in [0, 2, 5, 8, 10, 12, 15, 18, 20]:
    fwe = abs(fw[n] - exact[n]) / abs(exact[n])
    bwe = abs(bw[n] - exact[n]) / abs(exact[n])
    print(f"{n:>4} {exact[n]:>16.6e} {fw[n]:>18.6e} {fwe:>16.3e} "
          f"{bw[n]:>16.6e} {bwe:>16.3e}")

输出：

   n              精确值               前向递推           前向相对误差             后向递推           后向相对误差
   0     7.651977e-01       7.651977e-01        0.000e+00     7.651977e-01        1.451e-16
   2     1.149035e-01       1.149035e-01        0.000e+00     1.149035e-01        1.208e-16
   5     2.497577e-04       2.497577e-04        8.030e-12     2.497577e-04        1.085e-15
   8     9.422344e-08       9.422672e-08        3.480e-05     9.422344e-08        4.214e-16
  10     2.630615e-10       1.199871e-09        3.561e+00     2.630615e-10        3.931e-16
  12     4.999718e-13       4.101109e-07        8.203e+05     4.999718e-13        4.039e-15
  15     2.297532e-17       7.130539e-03        3.104e+14     2.297532e-17        1.744e-15
  18     5.880345e-22       2.320072e+02        3.945e+23     5.880345e-22        2.879e-15
  20     3.873503e-25       3.168944e+05        8.181e+29     3.873503e-25        1.897e-15

两列误差的走向完全相反。前向递推的相对误差按指数增长：$n=5$ 时 $8.03\times {10}^{-12}$（还能用），$n=8$ 时 $3.48\times {10}^{-5}$（已经损失五位有效数字），$n=10$ 时 $3.561$（误差超过真值本身），$n=20$ 时 $8.18\times {10}^{29}$（结果与真值毫无关系）。后向递推的相对误差在所有 $n$ 上都在 ${10}^{-15}$ 到 ${10}^{-16}$ 之间，达到双精度的极限。

失效的机制可以从递推的两个解看出来。$J_n(1)$ 随 $n$ 迅速衰减（从 $7.65\times {10}^{-1}$ 降到 $3.87\times {10}^{-25}$，跨越 $24$ 个数量级），而递推的另一个解 $Y_n(1)$ 随 $n$ 迅速增长；前向递推的放大因子 $\frac{2n}{x}$ 在 $n=20$ 时是 $40$，误差中的 $Y_n$ 分量被逐次放大，$n$ 步之后累积 ${40}^n$ 量级的放大。后向递推从 $n=60$ 向下进行，方向与 $Y_n$ 的增长方向相反，$Y_n$ 分量被逐次压缩，剩下的只是纯的 $J_n$ 分量，因此结果稳定。归一化用的是 $J_0+2\sum _{k\ge 1}{J}_{2k}=1$ 这条恒等式（对所有实数 $x$ 成立），它把「任意标度」的因素消掉。

实验四：奇异摄动与边界层

比较外解、复合展开与精确解，并对比三种有限差分格式在边界层问题上的表现。

import math
import numpy as np

print("eps*y'' + y' = x，y(0) = 0，y(1) = 2")
print("精确解 y = x^2/2 - eps*x + A - A e^(-x/eps)，A = (1.5 + eps)/(1 - e^(-1/eps))")
print("外解 y_out = x^2/2 + 1.5 + eps*(1 - x)，复合展开 y_comp = y_out - (1.5 + eps) e^(-x/eps)")

def exact_sol(xq, eps):
    e = math.exp(-1.0 / eps)
    A = (1.5 + eps) / (1.0 - e)
    return 0.5 * xq ** 2 - eps * xq + A - A * np.exp(-xq / eps)

def outer(xq, eps):
    return 0.5 * xq ** 2 + 1.5 + eps * (1 - xq)

def comp(xq, eps):
    return outer(xq, eps) - (1.5 + eps) * np.exp(-xq / eps)

xq = np.linspace(0, 1, 200001)
print(f"{'eps':>8} {'外解的最大误差':>18} {'复合展开的最大误差':>20}")
for eps in [0.1, 0.05, 0.02, 0.01]:
    ex = exact_sol(xq, eps)
    print(f"{eps:>8} {float(np.max(np.abs(outer(xq, eps) - ex))):>18.6f} "
          f"{float(np.max(np.abs(comp(xq, eps) - ex))):>20.3e}")

def solve_fd(eps, N, scheme):
    h = 1.0 / N
    A = np.zeros((N + 1, N + 1))
    b = np.zeros(N + 1)
    A[0, 0] = 1.0
    A[N, N] = 1.0
    b[N] = 2.0
    c = eps / h ** 2
    beta = math.exp(-h / eps)
    for i in range(1, N):
        xi = i * h
        if scheme == "central":
            A[i, i - 1] = c - 1.0 / (2 * h)
            A[i, i] = -2 * c
            A[i, i + 1] = c + 1.0 / (2 * h)
            b[i] = xi
        elif scheme == "upwind":
            A[i, i - 1] = c - 1.0 / h
            A[i, i] = -2 * c + 1.0 / h
            A[i, i + 1] = c
            b[i] = xi
        else:                                   # 指数拟合（Scharfetter-Gummel）
            A[i, i - 1] = c * beta
            A[i, i] = -c * (1 + beta)
            A[i, i + 1] = c
            b[i] = xi * (eps / h) * (1 - beta)
    return np.linspace(0, 1, N + 1), np.linalg.solve(A, b), beta

eps = 0.01
print(f"{'N':>6} {'h':>8} {'h/(2 eps)':>10} {'中心差分误差':>14} {'斜率变号次数':>14} "
      f"{'指数拟合误差':>14} {'迎风误差':>12}")
for N in [20, 50, 100, 500, 2000]:
    h = 1.0 / N
    xg, uc, _ = solve_fd(eps, N, "central")
    xg, ug, beta = solve_fd(eps, N, "sg")
    try:
        xg, uw, _ = solve_fd(eps, N, "upwind")
        uw_err = f"{float(np.max(np.abs(uw - exact_sol(xg, eps)))):>12.6f}"
    except np.linalg.LinAlgError:
        uw_err = f"{'矩阵奇异':>12}"
    ex = exact_sol(xg, eps)
    du = np.diff(uc)
    osc = int(np.sum(du[:-1] * du[1:] < 0))
    print(f"{N:>6} {h:>8.4f} {h / (2 * eps):>10.2f} "
          f"{float(np.max(np.abs(uc - ex))):>14.6f} {osc:>14} "
          f"{float(np.max(np.abs(ug - ex))):>14.6f} {uw_err}")

输出：

eps*y'' + y' = x，y(0) = 0，y(1) = 2
精确解 y = x^2/2 - eps*x + A - A e^(-x/eps)，A = (1.5 + eps)/(1 - e^(-1/eps))
外解 y_out = x^2/2 + 1.5 + eps*(1 - x)，复合展开 y_comp = y_out - (1.5 + eps) e^(-x/eps)
     eps            外解的最大误差            复合展开的最大误差
     0.1           1.600000            7.264e-05
    0.05           1.550000            3.195e-09
    0.02           1.520000            4.441e-16
    0.01           1.510000            4.441e-16
     N        h  h/(2 eps)         中心差分误差         斜率变号次数         指数拟合误差         迎风误差
    20   0.0500       2.50       0.657317              6       0.014469     1.857500
    50   0.0200       1.00       0.204356              0       0.002948         矩阵奇异
   100   0.0100       0.50       0.052165              0       0.000773     0.550548
   500   0.0020       0.10       0.001860              0       0.000031     0.060039
  2000   0.0005       0.03       0.000116              0       0.000002     0.014026

第一张表验证渐近展开的构造。外解的最大误差是 $1.600$、$1.550$、$1.520$、$1.510$，随 $\epsilon$ 减小趋向 $1.5$——误差不消失，位置固定在 $x=0$ 的边界层内（外解在 $x=0$ 处的值是 $1.5+\epsilon$，而真值是 $0$，差值约 $1.5$）。复合展开的最大误差是 $7.264\times {10}^{-5}$（$\epsilon =0.1$）、$3.195\times {10}^{-9}$（$\epsilon =0.05$）、$4.441\times {10}^{-16}$（$\epsilon =0.02$ 与 $0.01$，即机器精度）。$\epsilon$ 减半时误差降到约 $\frac{1}{2\times {10}^4}$，与「剩余误差是 $e^{-\text{常数}/\epsilon }$ 量级」的理论预期一致（$\epsilon$ 减半使 $e^{-1/\epsilon }$ 平方）。这条数据说明匹配渐近展开在边界层问题上的有效性：两个简单的解析表达式（外解与内解）就能把误差压到指数级小，代价是必须正确识别层的位置与宽度。

第二张表给出格式选择的三条结论。中心差分在 $h>2\epsilon$ 时出现振荡：$N=20$（$h/(2\epsilon )=2.5$）时最大误差 $0.657317$，解内部有 $6$ 处斜率变号；加密到 $N=50$（$h/(2\epsilon )=1.0$）后斜率变号消失，误差降到 $0.204356$，再加密到 $N=2000$ 时误差 $0.000116$。迎风格式在 $h=2\epsilon$（$N=50$）处矩阵奇异（对角元 $\frac{\epsilon }{h^2}-\frac{1}{h}=0$），求解失败；$h<2\epsilon$ 时格式可用但精度只有一阶（$N=2000$ 时 $0.014026$，比中心差分差两个数量级）。指数拟合格式在所有网格上误差最小（$N=20$ 时 $0.014469$，$N=2000$ 时 $0.000002$），斜率变号次数全部为零，且不依赖 $h<2\epsilon$ 这个条件。

三者的对照给出边界层问题的实践规则：当网格无法分辨层的宽度时，格式的单调性比形式精度更重要。中心差分的二阶精度只在 $h<2\epsilon$ 时兑现（$\epsilon ={10}^{-4}$ 时需要 $N>50000$）；指数拟合格式把方程在单元内的精确解的信息放进系数里，在 $h\gg \epsilon$ 时仍给出无振荡、量级正确的结果。工程上的常见组合是用指数拟合或迎风保证全局稳定，再在需要高精度的区域局部加密网格或改用高阶格式。

5.12.7 本节小结

渐近分析的三层内容

本节把渐近分析与特殊函数组织为三层。渐近方法层给出三件工具。Laplace 方法处理 $n$ 大时的指数型积分：主项 $e^{nf(x_0)}\sqrt{2\pi /(n|f^{″}(x_0)|)}$（内部最大值）或 $\frac{e^{nf(a)}}{n|f^{\prime }(a)|}$（端点最大值），修正项按 $1/n$ 的幂展开（实验一：$I(a)={\int }_{-\pi }^{\pi }{e}^{a\cos x}dx$ 的主项相对误差从 $a=1$ 的 $1.435\times {10}^{-1}$ 降到 $a=20$ 的 $6.395\times {10}^{-3}$，加 $1/(8a)$ 后按 $1/a^2$ 下降；$I(n)={\int }_0^1{e}^{-n{x}^2}dx$ 的主项误差按 $\frac{e^{-n}}{\sqrt{\pi n}}$ 衰减，最后一列从 $0.899$ 升到 $0.977$ 趋于 $1$）。Euler-Maclaurin 公式把求和化为积分加 Bernoulli 数修正：调和数 $H_N=\ln N+\gamma +\frac{1}{2N}-\frac{1}{12N^2}+\frac{1}{120N^4}-\frac{1}{252N^6}+\cdots$ 的四级误差实测为 $O(1/N)$、$O(1/N^2)$、$O(1/N^4)$、$O(1/N^6)$，系数逐项与 $\frac{1}{2}$、$\frac{1}{12}$、$\frac{1}{120}$、$\frac{1}{252}$ 相符。鞍点法把 Laplace 方法推广到复平面，与概率论的大偏差原理（Legendre 变换给出的速率函数）是同一个机制。特殊函数层给出三类函数的渐近形态（Stirling 公式、$\mathrm{erfc}$ 的尾部展开、Bessel 函数的 $x^{-1/2}$ 振幅振荡），并给出递推方向决定稳定性这条实用原则：Bessel 函数 $J_n(1)$ 的前向递推相对误差从 $n=5$ 的 $8.03\times {10}^{-12}$ 增长到 $n=20$ 的 $8.18\times {10}^{29}$（放大因子约 $(2n/x{)}^n$），后向递推（Miller 算法加归一化）在全部 $n$ 上保持 ${10}^{-15}$ 量级。摄动层给出匹配渐近展开的流程：外解（忽略最高阶导数项）、内解（放大坐标）、匹配条件、复合展开；实验四的数据：只用外解的最大误差是 $1.5$ 量级且不随 $\epsilon$ 消失（误差集中在边界层），复合展开的误差是 $7.264\times {10}^{-5}$（$\epsilon =0.1$）、$3.195\times {10}^{-9}$、$4.441\times {10}^{-16}$，随 $\epsilon$ 减小按指数下降；数值求解时中心差分在 $h>2\epsilon$ 处出现振荡（$N=20$ 时 $6$ 处斜率变号），迎风格式在 $h=2\epsilon$ 处矩阵奇异，指数拟合格式（Scharfetter-Gummel，系数恒正）在全部网格上保持单调且误差最小。

渐近分析给出了「参数极端时解的主项」，下一节作为选读把常微分方程与动力系统的视角补齐：解的存在唯一性、相平面与稳定性分析、数值方法的稳定性与刚性方程。这些内容与本章的其余部分关系密切：存在性用压缩映射原理（5.1 节），稳定性分析用线性化与 Jacobian 的谱（5.8 节），数值方法的稳定性用与本节相同的单调性与放大因子分析。

练习题

第 1 题 概念推导

(a) 用 Laplace 方法求 ${\int }_0^1{x}^n{e}^{-x}dx$ 的主项（提示：写成 ${\int }_0^1{e}^{n\ln x-x}dx$ 并注意最大值的位置）。(b) 写出 $\epsilon y^{″}+y^{\prime }=1$、$y(0)=0$、$y(1)=2$ 的外解与边界层位置，并说明为什么层出现在 $x=0$ 而不是 $x=1$。(c) 说明为什么 Bessel 函数的前向递推不稳定，而后向递推稳定。

参考答案

(a) 写成 $I_n={\int }_0^1{e}^{n\ln x-x}dx$，令 $f(x)=\ln x-\frac{x}{n}$（把 $e^{-x}$ 的部分归入 $f$）。$n$ 大时 $f^{\prime }(x)=\frac{1}{x}-\frac{1}{n}>0$ 在 $(0,n)$ 上成立，$f$ 在 $[0,1]$ 上单调递增，最大值在端点 $x=1$ 处取到，$f(1)=-\frac{1}{n}$，$f^{\prime }(1)=1-\frac{1}{n}\approx 1$。属于端点最大值情形：

$$I_n\approx \frac{e^{nf(1)}}{n|f^{\prime }(1)|}=\frac{e^{-1}}{n}\cdot \frac{1}{1-1/n}\approx \frac{e^{-1}}{n}$$

精确值可以用不完全 $\mathrm{\Gamma }$ 函数写出（$I_n=\gamma (n+1,1)$），主项 $\frac{1}{en}$ 与它的相对误差是 $O(1/n)$。这条推导说明端点情形的两个要点：最大值在端点时衰减是 $1/n$ 而不是 $n^{-1/2}$；端点的导数 $f^{\prime }(a)$ 出现在分母中，若 $f^{\prime }(a)=0$（端点也是平坦极值）则需要更高阶的处理。

(b) 忽略 $\epsilon y^{″}$ 得 $y_1^{\prime }=1$，即 $y_1=x+C$。用哪个边界定 $C$？取 $x=1$ 的边界条件 $y(1)=2$，得 $C=1$，外解是 $y_1=x+1$。层出现在 $x=0$，判断依据是：一阶方程的解只满足一个边界条件，另一个边界条件对应的位置就是层的所在。具体地，$y_1(0)=1\ne y(0)=0$，差值 $1$ 必须在 $x=0$ 附近的薄层内被消掉，层宽由 $\xi =\frac{x}{\epsilon }$ 定出为 $O(\epsilon )$。内解是 $Y(\xi )=A+B{e}^{-\xi }$，配 $Y(0)=0$ 与匹配条件 $Y(\mathrm{\infty })=y_1(0)=1$ 得 $A=1$、$B=-1$，复合展开 $y_{\text{comp}}=x+1-e^{-x/\epsilon }$。若把两个边界条件反过来（$y(0)=2$、$y(1)=0$），层会移到 $x=1$，因为此时外解 $y_1=x-1$ 在 $x=1$ 处给出 $0$ 而边界条件是 $y(1)=0$ 刚好吻合，矛盾转到 $x=0$。判断规则：外解只能满足一个边界条件，被放弃的那个边界条件的位置就是层的位置；对流项系数为正（$y^{\prime }$ 的系数为 $+1$）时层在下游端（$x=1$）还是上游端取决于外解的选取，但通常层出现在「信息流入的边界」——本例中 $y^{\prime }$ 的系数为正意味着信息从左向右传播，因此层出现在 $x=0$（上游端）。

(c) 递推 $J_{n+1}=\frac{2n}{x}{J}_n-J_{n-1}$ 有两个线性无关的解：$J_n(x)$（随 $n$ 迅速衰减）与 $Y_n(x)$（随 $n$ 迅速增长）。任何数值计算都带有舍入误差，它相当于在解的表示中加入了一个小量的 $Y_n$ 分量（因为递推的通解是两者的线性组合）。前向递推时，$Y_n$ 分量每步被乘以约 $\frac{2n}{x}$ 的因子放大，$n$ 步之后放大 $(2n/x{)}^n$ 倍，从 ${10}^{-16}$ 的初始误差增长到 $O(1)$ 只需要几步（$x=1$ 时 $n=10$ 处的放大因子已经远大于 ${10}^{16}$）。后向递推时，$Y_n$ 分量被同样因子缩小，几步之后它相对 $J_n$ 可以忽略，剩下的几乎纯是 $J_n$ 分量。

代价是后向递推的起始值是任意的（不知道 $J_N$ 的准确值），但递推是线性的，任意起点给出的都是「$J_n$ 与 $Y_n$ 的某个组合」；由于 $Y_n$ 分量被压缩掉，结果与 $J_n$ 只差一个常数因子，用一条已知的恒等式（如 $J_0+2\sum J_{2k}=1$）归一化即可。这条方法称 Miller 算法，适用于所有「希望得到衰减解」的递推。判断方向的通用方法：比较递推系数的大小，若 $|a_n|>1$（这里 $a_n=\frac{2n}{x}$），则前向放大、后向缩小；再确认希望得到的解是衰减的还是增长的，选与它一致的递推方向。

第 2 题 计算推理

(a) 用 Euler-Maclaurin 公式估计 $\sum _{k=1}^N\sqrt{k}$ 的主项与第一修正项。(b) 用 Laplace 方法求 $(\genfrac{}{}{0ex}{}{2n}{n})$ 的主项。(c) 估计 ${\int }_0^{\mathrm{\infty }}\frac{e^{-x}}{1+x/n}dx$ 在 $n$ 大时的主项与第一修正项。

参考答案

(a) 取 $f(x)=\sqrt{x}$，Euler-Maclaurin 公式给出

$$\sum _{k=1}^N\sqrt{k}={\int }_1^N\sqrt{x}dx+\frac{\sqrt{1}+\sqrt{N}}{2}+\frac{B_2}{2!}(f^{\prime }(N)-f^{\prime }(1))+\cdots$$

计算各项：${\int }_1^N\sqrt{x}dx=\frac{2}{3}(N^{3/2}-1)$；$\frac{\sqrt{1}+\sqrt{N}}{2}=\frac{1}{2}+\frac{\sqrt{N}}{2}$；$f^{\prime }(x)=\frac{1}{2\sqrt{x}}$，$B_2=\frac{1}{6}$，故 $\frac{1}{12}(\frac{1}{2\sqrt{N}}-\frac{1}{2})$。合并：

$$\sum _{k=1}^N\sqrt{k}=\frac{2}{3}{N}^{3/2}+\frac{1}{2}{N}^{1/2}+\zeta (-\frac{1}{2})+\frac{1}{24}{N}^{-1/2}+O\left(N^{-3/2}\right)$$

（常数项 $\zeta (-\frac{1}{2})\approx -0.2079$ 由把 $N\to \mathrm{\infty }$ 时的极限整理得到，也可以通过直接数值求和核对）。主项 $\frac{2}{3}{N}^{3/2}$ 与第一修正项 $\frac{1}{2}{N}^{1/2}$ 的量级比是 $O(N^{-1})$，$N={10}^4$ 时修正占主项的 $0.75\mathrm{\%}$。这条展开在算法分析中用于估计「分块求和的代价」这类量。

(b) $(\genfrac{}{}{0ex}{}{2n}{n})=\frac{(2n)!}{(n!{)}^2}$，用 Stirling 公式 $\ln m!=m\ln m-m+\frac{1}{2}\ln (2\pi m)+O(1/m)$：

$$\ln (\genfrac{}{}{0ex}{}{2n}{n})=\ln (2n)!-2\ln n!=2n\ln (2n)-2n+\frac{1}{2}\ln (4\pi n)-2(n\ln n-n+\frac{1}{2}\ln (2\pi n))+O(1/n)$$

简化：$2n\ln (2n)-2n\ln n=2n\ln 2$，$-2n+2n=0$，$\frac{1}{2}\ln (4\pi n)-\ln (2\pi n)=-\frac{1}{2}\ln (\pi n)$。因此

$$\ln (\genfrac{}{}{0ex}{}{2n}{n})=2n\ln 2-\frac{1}{2}\ln (\pi n)+O(1/n)⟹(\genfrac{}{}{0ex}{}{2n}{n})\sim \frac{4^n}{\sqrt{\pi n}}$$

数值核对：$n=10$ 时 $\frac{4^{10}}{\sqrt{10\pi }}=\frac{1048576}{5.605}=187078$，而 $(\genfrac{}{}{0ex}{}{20}{10})=184756$，相对误差 $1.3\mathrm{\%}$；$n=100$ 时误差降到 ${10}^{-3}$ 量级。这条结果在信息论中用于估计「长度为 $n$ 的二进制串中熵最大的那类」的规模（第四章的组合数估计），也是熵 $H=1$ 比特的体现：$4^n=2^{2n}$ 恰好是 $2n$ 个比特的全部可能数。

(c) 把被积函数写成 $e^{-x}{e}^{-\ln (1+x/n)}$ 并展开对数：$\ln (1+x/n)=\frac{x}{n}-\frac{x^2}{2n^2}+\cdots$，因此

$${\int }_0^{\mathrm{\infty }}\frac{e^{-x}}{1+x/n}dx={\int }_0^{\mathrm{\infty }}{e}^{-x(1+\frac{1}{n})+\frac{x^2}{2n^2}-\cdots }dx$$

主项取到 $O(1/n)$：${\int }_0^{\mathrm{\infty }}{e}^{-x(1+1/n)}dx=\frac{1}{1+1/n}=1-\frac{1}{n}+\frac{1}{n^2}-\cdots$。因此主项是 $1$、第一修正项是 $-\frac{1}{n}$。这也是 ${\int }_0^{\mathrm{\infty }}\frac{e^{-x}}{1+x/n}dx$ 的部分分式展开（$=n{e}^n{E}_1(n)$，$E_1$ 是指数积分）的渐近形式，与 $\frac{1}{n}$ 的幂次展开一致。这条计算的实用含义：把分母展开成 $1+\frac{x}{n}$ 的形式时，修正项的量级由 $\frac{E[X]}{n}$ 决定（$X$ 服从指数分布，$E[X]=1$），因此 $n$ 必须远大于被积函数的一阶矩才能保证主项有效。

第 3 题 代码验证

用本节代码做三组实验并解释：(a) 在实验一里把积分换成 $I(n)={\int }_0^1{e}^{-n{x}^4}dx$，观察主项从 $n^{-1/2}$ 变成什么；(b) 在实验三里把 $x$ 从 $1$ 改成 $10$，观察前向递推的失稳点移动到哪个 $n$；(c) 在实验四里把 $\epsilon$ 改成 $0.05$，重新确定中心差分保持无振荡的最大步长，并与 $2\epsilon$ 的预测比较。

参考答案

(a) 作换元 $t=n{x}^4$：$x=(t/n{)}^{1/4}$，$dx=\frac{1}{4}{t}^{1/4-1}{n}^{-1/4}dt$，因此

$$I(n)={\int }_0^1{e}^{-n{x}^4}dx=\frac{1}{4}{n}^{-1/4}{\int }_0^n{t}^{-3/4}{e}^{-t}dt\to \frac{1}{4}\mathrm{\Gamma }\left(\frac{1}{4}\right)n^{-1/4}$$

主项从 $n^{-1/2}$ 变成 $n^{-1/4}$，即 $O(n^{-1/(2m)})$ 中的 $m=2$（$f(x)=-x^4$ 在起点的展开始项是 $x^4$）。数值核对：$\frac{1}{4}\mathrm{\Gamma }(0.25)=\frac{3.6256}{4}=0.9064$，$n=100$ 时 $0.9064/{100}^{0.25}=0.9064/3.1623=0.2866$，而精确积分的量级与此相符。这条规律的实用含义是：主项的幂次由最大值附近的平坦程度决定，知道 $f$ 的展开始项就能预测幂次，不必逐次计算。

(b) 递推的放大因子是 $\frac{2n}{x}$，初始误差约 ${10}^{-16}$，失稳条件是 ${\left(\frac{2n}{x}\right)}^n$ 超过 ${10}^{16}$。$x=10$ 时放大因子是 $\frac{2n}{10}=\frac{n}{5}$，比 $x=1$ 时小 $10$ 倍，因此需要更大的 $n$ 才能积累到同样的误差。粗略估计：$x=1$ 时 $n=10$ 处已经失控（放大因子 $20$，${20}^{10}={10}^{13}$，与 ${10}^{16}$ 同量级）；$x=10$ 时要求 ${\left(\frac{n}{5}\right)}^n\sim {10}^{16}$，取 $n=16$ 得 ${3.2}^{16}=1.8\times {10}^8$，$n=25$ 得 $5^{25}=3\times {10}^{17}$，因此失稳点大致后移到 $n=22$ 到 $25$。结论：$x$ 越大（相对 $n$ 而言）递推越稳定，因为 $J_n(x)$ 在 $n<x$ 的区间内并不衰减（振荡区），前向递推在那里是稳定的；失稳只发生在 $n$ 超过 $x$ 之后的衰减区。

(c) $\epsilon =0.05$ 时中心差分的稳定条件 $h<2\epsilon =0.1$ 对应 $N>10$。取 $N=10$ 是临界（$h=0.1=2\epsilon$），所以无振荡的最小网格数是 $N=12$ 或 $16$（$h=0.083$ 或 $0.0625$）。运行实验代码可以核对：$N=20$（$h=0.05<0.1$）时斜率变号次数应当为 $0$，这与 $\epsilon =0.01$ 时 $N=20$ 出现 $6$ 次变号形成对照。这条练习说明网格要求随 $\epsilon$ 线性变化：$\epsilon$ 减半，保持无振荡所需的网格数加倍，因此 $\epsilon ={10}^{-6}$ 时需要 $N$ 在 ${10}^7$ 量级——这正是边界层问题必须使用自适应网格或指数拟合格式的量化理由。

第 4 题 综合应用

某数值程序需要计算 $I(n)={\int }_0^1\frac{e^{-nx}}{1+x}dx$（$n$ 从 $1$ 到 ${10}^6$ 不等，要求相对误差低于 ${10}^{-6}$），并需要给出每个 $n$ 处结果的误差估计。请设计计算方案，比较可选的路线（数值积分、Laplace 方法、级数展开），说明各自的误差来源与适用范围。

参考答案

分析积分结构：$n$ 小时被积函数在整个区间上平缓，$n$ 大时被积函数在 $x=0$ 附近急剧下降（衰减尺度 $1/n$）。因此单一方法不能覆盖全部参数范围，需要分段处理。分界点由「衰减尺度与区间长度之比」决定：$n\lesssim 10$ 时衰减尺度 $\frac{1}{n}$ 与 $1$ 可比，$n\gtrsim 10$ 时衰减尺度远小于区间。

路线一：数值积分（自适应 Gauss）。用自适应求积（例如 Gauss-Kronrod 的 G7-K15 对，误差由两个结果的差估计）在 $[0,1]$ 上计算。误差来源是求积误差，可用自适应细分控制：对 $n\lesssim {10}^3$ 的情形，被积函数光滑，自适应求积用几百次求值就能达到 ${10}^{-10}$ 的精度，误差估计直接来自 G7 与 K15 的差。局限：$n$ 很大时需要把区间细分到 $O(1/n)$ 量级才能分辨初始的急剧下降，求值次数随 $n$ 增长（约 $O(\log n)$ 到 $O(n^{1/2})$ 量级），$n={10}^6$ 时可能不划算。

路线二：级数展开（渐近级数）。把 $\frac{1}{1+x}$ 展开成 $1-x+x^2-\cdots$ 并逐项积分：

$$I(n)=\sum _{k=0}^{\mathrm{\infty }}(-1{)}^k{\int }_0^1{x}^k{e}^{-nx}dx$$

而 ${\int }_0^1{x}^k{e}^{-nx}dx=\frac{k!}{n^{k+1}}(1-e^{-n}\sum _{j=0}^k\frac{n^j}{j!})$。$n$ 大时主项是 $\frac{1}{n}$（$k=0$ 项），第一修正项是 $-\frac{1}{n^2}$（$k=1$ 项），因此

$$I(n)=\frac{1}{n}-\frac{1}{n^2}+\frac{2}{n^3}-\frac{6}{n^4}+\cdots =\frac{1}{n}(1-\frac{1}{n}+\frac{2}{n^2}-\frac{6}{n^3}+\cdots )$$

截到 $k=K$ 时相对误差是 $\frac{(K+1)!}{n^{K+1}}$ 量级。$n={10}^6$ 时取两项（$K=1$）即得 $2\times {10}^{-6}$ 量级，取三项（$K=2$）后为 $6\times {10}^{-12}$，远超要求。$n=10$ 时情形完全不同：项 $(K+1)!/{10}^{K+1}$ 随 $K$ 先减后增（比值是 $\frac{K+2}{10}$，$K>8$ 之后开始增大），最优截断处的相对误差是 $9!/{10}^9=3.6\times {10}^{-4}$，达不到 ${10}^{-6}$ 的要求。这正是渐近级数的特征：它在大参数处极其精确，在中等参数处精度有上限（由最小项决定），继续加项还会变差。因此级数展开只适用于 $n\ge {10}^3$ 的情形，且必须「截到最小项为止」。

路线三：换元 + 数值积分。把积分变量改成 $t=nx$：$I=\frac{1}{n}{\int }_0^n{e}^{-t}\frac{dt}{1+t/n}$。$n$ 大时积分上限的贡献可以忽略，被积函数在 $t\in [0,\mathrm{\infty })$ 上光滑且指数衰减，自适应求积非常高效（几十次求值即可），且误差估计仍然来自求积公式。这条路线结合了数值的稳健与解析的形式，适用于全部参数范围：$n$ 小时用原区间加自适应细分，$n$ 大时用换元后的形式。

推荐方案：分段实施。$n\le {10}^2$ 时用自适应 Gauss 求积（原区间或换元后都可，换元后更稳）；$n>{10}^2$ 时用换元后的形式加自适应求积，同时在 $n>{10}^3$ 时用级数展开（取到最小项）作为独立核对，两条路线的差作为误差估计。误差估计的设计：求积路线用 G7-K15 的差（可靠且方便），级数路线用截断项的量级（$\frac{(K+1)!}{n^{K+1}}$），两者在重叠区（$n$ 在 ${10}^2$ 到 ${10}^3$）比较，确认一致后再报告。报告中给出每个 $n$ 所用的方法与误差估计，这对「$n$ 跨越六个数量级」的批量计算是必需的——没有任何一种方法对全部参数范围都最优，方法的切换点与切换处的交叉验证是方案设计的一部分。

常见错误

错误 1 · 在参数不够大时使用渐近公式

原因：渐近展开的每一项都假设「展开参数是小量」，参数不够大时误差很大，甚至主项的符号或量级都可能错。本节实验一第二组的数据给出量化对照：$a=1$ 时主项 ${\int }_{-\pi }^{\pi }{e}^{a\cos x}dx$ 的相对误差是 $1.435\times {10}^{-1}$（$a=20$ 时降到 $6.395\times {10}^{-3}$）；即使加上 $1/(8a)$ 修正，$a=1$ 时的相对误差仍是 $3.639\times {10}^{-2}$。$a=1$ 时修正项 $1/(8a)=0.125$ 已经不是小量，展开失去意义。另一类问题是混淆不同的渐近区域：${\int }_0^1{e}^{-n{x}^2}dx$ 的主项 $\frac{\sqrt{\pi }}{2\sqrt{n}}$ 需要 $n$ 大，而 ${\int }_0^1{e}^{-n{x}^4}dx$ 的主项是 $\frac{\mathrm{\Gamma }(1/4)}{4}{n}^{-1/4}$，两者的幂次不同，用错公式会得到完全不同的量级。

解决：使用渐近公式前检查参数是否落在适用区间（一般要求展开参数 $\le 0.1$ 才能保证 $1\mathrm{\%}$ 的精度），并给出误差估计（用下一项的量级，或者与数值积分对照）。报告中写明「渐近公式 + 适用范围」，对小参数的情形改用数值方法或精确公式。若参数跨越多个量级，分段使用不同方法并给出切换点的验证（如练习 4 的方案）。

错误 2 · 把发散的渐近级数当作收敛级数求和

原因：渐近级数的系数（Bernoulli 数、阶乘项）增长很快，级数对固定的参数最终发散。本节实验二的数据给出了边界现象：$N=1000$ 之后第三、四列修正的误差落到 ${10}^{-15}$ 并保持不变，这是舍入误差的下限；若继续加入更多项（如 $\frac{1}{252N^6}$ 之后的项），误差会先减小后回升，最终发散（对调和数，Bernoulli 数的增长使级数在项数约 $\pi N$ 处开始发散）。把渐近级数当成 Taylor 级数那样「加得越多越准」会得到错误的结论，甚至出现「误差随项数增加而上升」的结果。

解决：遵循「截到最小项为止」的原则：逐项累加并在误差开始回升时停止（与 5.4 节数值微分的最优步长是同一类权衡）。用余项的估计给出误差界（对交错型的渐近级数，误差小于第一个被舍弃项）。对需要高精度的场合，改用收敛的展开（如 Taylor 级数、连分式、Padé 逼近）或者数值方法。报告中说明所用的项数与误差估计，不要报告「加到第 $k$ 项」而不说明为什么停在那里。

错误 3 · 选择与解的衰减方向相反的递推方向

原因：三项递推有两个解，一个增长一个衰减，递推方向决定了误差的放大或压缩。本节实验三的数据给出极端例子：Bessel 函数 $J_n(1)$ 的前向递推在 $n=20$ 时的相对误差是 $8.18\times {10}^{29}$（结果完全没有意义），而后向递推保持在 $1.9\times {10}^{-15}$。同类问题出现在很多递推计算中：连带 Legendre 函数、球谐函数、Airy 函数、Chebyshev 展开的系数递推、三对角方程组的部分算法都可能有这个问题。隐藏性在于公式本身完全正确，错误只体现在数值结果上，而且往往「前几步看起来正常」，问题在后面才暴露。

解决：先判断递推的两个解中哪个是你需要的、哪个是增长的（比较递推系数的量级）；若要得到衰减解，用后向递推加归一化（Miller 算法）；若两者都不衰减（临界情形），需要更细致的分析（例如用 Chebyshev 型的方法或改写成稳定的二阶形式）。实现时做「方向对照实验」：用两个方向各算一次，看结果是否一致，不一致说明方向选错。对已经算出的结果，用已知恒等式或特殊点上的已知值核对。

错误 4 · 在边界层问题中使用常规差分格式而不检查网格条件

原因：含小参数的微分方程在边界层内的解变化尺度是 $O(\epsilon )$，常规的中心差分要求 $h<2\epsilon$ 才能保持单调，否则产生非物理振荡。本节实验四的数据给出完整形态：$\epsilon =0.01$、$N=20$（$h/(2\epsilon )=2.5$）时中心差分的最大误差是 $0.657317$，解中有 $6$ 处斜率变号；迎风格式在 $h=2\epsilon$ 处矩阵奇异（对角元为零），直接求解失败。这两类现象（振荡与奇异）在调试时容易误判为「数据错误」或「程序 bug」，实际上是格式与网格不匹配。

解决：先估计解的变化尺度（由方程的小参数与系数决定，通常是 $\epsilon$ 或 $\epsilon$ 的幂），据此检查网格条件（对流扩散方程要求 $h<2\epsilon$，更高阶的方程有相应的条件）。不满足条件时按优先级选择对策：改用单调格式（指数拟合/Scharfetter-Gummel、迎风）、使用边界层拟合的网格（在层内几何加密）、或者做变量伸缩把层拉平。实现后用「解的单调性」与「网格加密时的收敛率」做检验：斜率变号次数应当为零，收敛率应当与格式的阶数相符。报告中给出 $\epsilon$、$h$、网格 Péclet 数与所用格式，这些信息决定了结果是否可信。