企业官网建设流程全解析

1. 一次拟合实验暴露的问题：普通多项式基在第 10 阶就崩了

先说个我印象很深的小事故。早年做一组球面测量数据的曲面拟合，图省事直接拿 $1, x, x^2, \dots, x^n$ 当基函数，构造最小二乘的正规方程求解系数。前 6 阶看起来都挺正常，拟合残差一路下降；阶数加到 10 之后，残差不但没降，系数反而开始正负跳动，量级从 $10^{-2}$ 一路飙到 $10^{5}$，好不容易降到 $10^{-4}$ 的误差又弹回去了。当时第一反应是数据有噪声，换了几批数据都一样，才意识到问题出在基函数本身。

这就是勒让德多项式登场的地方。它是定义在 $[-1, 1]$ 上的一族正交多项式，正交性让最小二乘的正规方程从"稠密的病态矩阵"退化成"对角矩阵"，每个系数的求解彼此独立，数值稳定性直接拉开一个数量级。它的来源有两个看似不相关的方向：一边是电动力学里球坐标系下解拉普拉斯方程的分离变量，另一边是函数逼近论里的最佳平方逼近。同一个多项式在两处出现，本质是同一个数学结构在不同问题里的投影。

这篇文章想做的事情很具体：把勒让德多项式从"课本上一个带积分符号的定义"变成你手上能直接调、能直接验、能直接判断"这里该用几阶"的工具。适合正在做数值拟合、球面/轴对称场计算、有限元或谱方法、以及需要做高精度数值积分的读者。数学基础不需要太深，会求导、会一点线性代数、能跑 Python 就够了，关键的推导我会把每一步的"为什么这么走"说清楚。

2. 单项式基的病态到底来自哪儿

2.1 从正规方程的矩阵结构看问题

最小二乘要解的是 $A^{\mathsf T}A \mathbf{c} = A^{\mathsf T}\mathbf{y}$，其中 $A$ 的列是各个基函数在采样点上的取值。如果基函数取 $1, x, x^2, \dots$，且在 $[-1,1]$ 上近似均匀采样，那么 $A^{\mathsf T}A$ 的元素近似正比于 $\int_{-1}^{1} x^{i+j} dx$，也就是 $\frac{2}{i+j+1}$ 乘以一个常数。这个矩阵就是希尔伯特矩阵的变体：它的条件数随阶数指数级恶化，$n=10$ 时已经在 $10^{13}$ 量级。

条件数是什么概念？它意味着输入数据的第 13 位有效数字被放大到影响解的第一位。双精度浮点只有约 16 位有效数字，等于说在 10 阶附近，你的系数解里基本只剩噪声。此时求解出来的系数在数学上"存在"，在数值上"没用"，振荡就是这么来的。

2.2 正交基带来的矩阵对角化

换一组基 $\phi_0, \phi_1, \dots$，如果它们在区间上满足 $\int \phi_i \phi_j w,dx = 0\ (i \neq j)$，那么 $A^{\mathsf T}A$ 的非对角元全部接近零，只剩对角元素 $h_i = \int \phi_i^2 w ,dx$。正规方程立刻退化成 $n+1$ 个独立的一元一次方程：

$$c_i = \frac{\int f(x)\phi_i(x)w(x),dx}{\int \phi_i^2(x) w(x),dx}$$

每个系数单独算，互不影响，也不存在矩阵求逆。条件数从 $10^{13}$ 变成 1，这就是正交基的全部意义。勒让德多项式正是取权重 $w(x) \equiv 1$、区间 $[-1,1]$ 时的那一族。

2.3 它在一堆正交多项式里站的位置

常见的正交多项式族其实不少，选哪个取决于区间和权重。我按实际使用频率整理了一张对照：

关键判断依据是问题的自然区间和自然权重。如果你的量本来就定义在 $[0,\infty)$ 且按 $e^{-x}$ 衰减，硬套勒让德多项式就得先做区间截断，截断误差会反过来吃掉正交性带来的精度收益。反过来，球坐标里的极角 $\theta$ 换元成 $x = \cos\theta$ 之后，积分自然落在 $[-1,1]$ 且权重恰好是 1，这时候勒让德就是唯一合适的选择。

3. 三种定义各管什么用：求导公式、递推与生成函数

3.1 罗德里格斯公式：理论推导的起点

教科书上最常见的形式是罗德里格斯公式：

$$P_n(x) = \frac{1}{2^n n!}\frac{d^n}{dx^n}\left[(x^2-1)^n\right]$$

这个式子在做理论推导时非常好用——要证明正交性、证明 $P_n$ 的 $n$ 个零点都在 $(-1,1)$ 内、推导递推关系，基本都从这里出发。思路是对方程 $(x^2-1)^n$ 反复求导 $n$ 次，得到一个 $n$ 次多项式，再验证它满足勒让德方程 $(1-x^2)y'' - 2xy' + n(n+1)y = 0$。

但要提醒一句：这个公式在代码里千万别直接用。$n=20$ 时 $2^{20} \cdot 20!$ 已经是天文数字，展开 $(x^2-1)^n$ 会产生严重的系数抵消，双精度下 $n$ 超过 15 左右就开始出现肉眼可见的误差。它是给人看的，不是给机器算的。

3.2 三项递推：工程实现唯一推荐路径

真正在代码里该用的是这一条：

$$(n+1)P_{n+1}(x) = (2n+1)x P_n(x) - n P_{n-1}(x)$$

从 $P_0(x)=1$、$P_1(x)=x$ 出发，一路推到任意阶，每步只用到一次乘法和一次减法，没有大数相除，没有系数抵消，数值上非常稳。前几项手算出来长这样：

$$P_0 = 1,\quad P_1 = x,\quad P_2 = \frac{3x^2-1}{2},\quad P_3 = \frac{5x^3-3x}{2},\quad P_4 = \frac{35x^4-30x^2+3}{8}$$

实测在 $[-1,1]$ 上，用递推算到 $n=100$ 时 $|P_n(x)|$ 仍在 $O(1)$ 量级（理论上 $|P_n(x)| \le 1$），说明没有发生数值爆炸。这也解释了为什么所有成熟数值库的底层实现都是递推而非求导公式。

3.3 生成函数与 $1/|\mathbf{r}-\mathbf{r}'|$ 的血缘关系

第三个定义是生成函数，它把整族多项式打包进一个解析表达式：

$$\frac{1}{\sqrt{1 - 2xt + t^2}} = \sum_{n=0}^{\infty} P_n(x),t^n,\qquad |t|<1$$

这个式子左边在物理里出现的频率高得惊人——它就是两点距离的倒数。设 $\mathbf{r}$ 与 $\mathbf{r}'$ 的夹角为 $\gamma$，长度分别为 $r$ 和 $r'$，那么当 $r' < r$ 时：

$$\frac{1}{|\mathbf{r}-\mathbf{r}'|} = \frac{1}{r}\sum_{n=0}^{\infty}\left(\frac{r'}{r}\right)^n P_n(\cos\gamma)$$

电动力学里的多极展开、静电场势能展开、甚至天体力学里的引力势摄动项，都是从这个式子长出来的。把 $x = \cos\gamma$、$t = r'/r$ 代进生成函数直接得到，不需要额外推导。这也是我建议记住生成函数的原因：它把"多项式"和"物理里的距离倒数"这两件事缝在了一起，看到 $1/|\mathbf{r}-\mathbf{r}'|$ 就该条件反射想到勒让德展开。

4. 正交区间与权重：为什么偏偏是 $[-1,1]$ 和常数 1

4.1 内积的定义与归一化常数

勒让德多项式的正交关系写成：

$$\int_{-1}^{1} P_m(x)P_n(x),dx = \frac{2}{2n+1}\delta_{mn}$$

注意右边不是 1，而是 $\frac{2}{2n+1}$。这个细节在推导最佳逼近系数时极其关键——很多人第一次写代码时忘了除以它，结果算出来的系数整体偏大，误差图上表现为"拟合曲线完全跟住了趋势但幅度不对"，看着像缩放错误，实际是漏归一化。

如果要做成标准正交（即范数为 1），需要再乘上 $\sqrt{\frac{2n+1}{2}}$。工程上我更倾向保持原始归一化，因为系数公式里那个 $\frac{2n+1}{2}$ 因子和正交关系里的分母正好配对，写起来更直观。真要标准化，就在输出层面统一处理。

4.2 区间变换的标准写法

实际问题很少恰好定义在 $[-1,1]$ 上。设原始区间是 $[a,b]$，标准变换是：

$$x = \frac{2t - (a+b)}{b-a},\qquad t \in [a,b]$$

反过来 $t = \frac{(b-a)x + (a+b)}{2}$。这里有个极易踩的坑：积分里的 $dx$ 也要跟着换。原积分 $\int_a^b f(t)P_n(x(t)),dt$ 换元后变成 $\frac{b-a}{2}\int_{-1}^{1} f(t(x))P_n(x),dx$，雅可比因子 $\frac{b-a}{2}$ 一旦漏掉，所有系数都会差一个固定倍数。这个错误隐蔽性很强，因为曲线形状完全正确，只是幅值不对，很容易被误判成"数据本身有偏差"。

我一般会在代码里把这个雅可比单独写成一个变量jac = (b - a) / 2，并且在计算完系数后做一次重构校验，直接对比 $f$ 和重构结果的最大偏差。校验不过就别往下走了。

4.3 奇偶性、零点分布与端点行为

掌握几个结构性质，能帮你在结果异常时快速定位问题：

• 奇偶性：$P_n(-x) = (-1)^n P_n(x)$。偶数阶是偶函数，奇数阶是奇函数。如果你的目标函数 $f$ 是偶函数，理论上的奇次项系数应当严格为零，实际算出来如果是 $10^{-16}$ 量级那就是浮点残差，如果到了 $10^{-3}$ 量级，说明区间变换或权重搞错了。
• 零点分布：$P_n$ 在 $(-1,1)$ 内有 $n$ 个互不相同的实零点，且关于原点对称。这些零点就是高斯-勒让德积分节点，后面还会用到。
• 端点值：$P_n(1) = 1$，$P_n(-1) = (-1)^n$。所有阶在端点都取到 $\pm 1$，这是级数在端点附近收敛慢的根源。

最后一条值得展开说。$|P_n(x)| \le 1$ 在整个区间恒成立，意味着每一项对结果的贡献上限就是系数本身。所以系数的衰减速度直接决定了需要多少项。如果 $f$ 是解析函数（比如 $e^x$），系数按阶乘级衰减，十来项就够；如果 $f$ 只有有限阶导数（比如 $|x|$），衰减只能是多项式级，需要几百项才能压住误差。

5. 球坐标下的轴对称边值问题：分离变量怎么落到它头上

5.1 拉普拉斯方程分离出勒让德方程

球坐标里，轴对称（与方位角 $\varphi$ 无关）的拉普拉斯方程是：

$$\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial V}{\partial r}\right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial V}{\partial\theta}\right) = 0$$

设 $V = R(r)\Theta(\theta)$ 代入，两边同乘 $r^2/(R\Theta)$，可以得到与角度有关的独立方程。做换元 $x = \cos\theta$，$\Theta(\theta)$ 变成 $y(x)$，方程化为：

$$(1-x^2)y'' - 2xy' + l(l+1)y = 0$$

这就是勒让德方程。它在 $x = \pm 1$（对应 $\theta = 0, \pi$ 两个极点）处有正则奇点，要求解在极点上有限，正好把 $l$ 约束成非负整数，对应的解就是 $P_l(x)$。第二类解 $Q_l$ 在端点发散，物理上被排除——这个排除不是人为选择，是物理边界条件逼出来的。

于是轴对称解的一般形式是：

$$V(r,\theta) = \sum_{l=0}^{\infty}\left(A_l r^l + B_l r^{-(l+1)}\right)P_l(\cos\theta)$$

$r^l$ 项在原点处有限，适合球内区域；$r^{-(l+1)}$ 项在无穷远处衰减，适合球外区域。选哪一组取决于你要解的区域。

5.2 边界条件匹配靠的就是正交性

假设球面 $r=a$ 上给定电位 $V_0(\theta)$，求球内电位。球内不能含 $r^{-(l+1)}$ 项（否则原点是奇点），所以：

$$V(r,\theta) = \sum_{l=0}^{\infty} A_l r^l P_l(\cos\theta)$$

在 $r=a$ 上要求 $\sum A_l a^l P_l(\cos\theta) = V_0(\theta)$。两边同乘 $P_m(\cos\theta)\sin\theta$ 并对 $\theta$ 从 $0$ 到 $\pi$ 积分，左边只有 $l=m$ 项存活，右边得到一个积分。整理后：

$$A_l a^l = \frac{2l+1}{2}\int_{-1}^{1} V_0(\arccos x),P_l(x),dx$$

这个式子是全套流程的核心。它之所以成立，全靠 $\int_{-1}^{1}P_l P_m,dx = \frac{2}{2l+1}\delta_{lm}$ 这一个性质——没有正交性，这个无穷级数的系数根本没法分离出来。

5.3 一个完整算例：上下半球电位相反的球壳

取具体的例子好说明问题。半径 $a$ 的球壳，上半球（$\theta < \pi/2$）电位 $+V_0$，下半球电位 $-V_0$。换成 $x = \cos\theta$，边界条件就是 $V_0(x) = V_0 \cdot \mathrm{sign}(x)$，一个奇函数。

偶次项系数应为零，只算奇次项：

$$A_l a^l = \frac{2l+1}{2}\cdot 2V_0\int_0^1 P_l(x),dx$$

手算前几项：$\int_0^1 P_1,dx = \frac{1}{2}$，$\int_0^1 P_3,dx = -\frac{1}{8}$，$\int_0^1 P_5,dx = \frac{1}{16}$。代进去得到：

$$A_1 a = \frac{3}{2}V_0,\qquad A_3 a^3 = -\frac{7}{8}V_0,\qquad A_5 a^5 = \frac{11}{16}V_0$$

所以球内电位是：

$$V(r,\theta) = V_0\left[\frac{3}{2}\left(\frac{r}{a}\right)P_1(\cos\theta) - \frac{7}{8}\left(\frac{r}{a}\right)^3 P_3(\cos\theta) + \frac{11}{16}\left(\frac{r}{a}\right)^5P_5(\cos\theta) - \cdots\right]$$

有个现成的检验方法：在 $\theta = 0$（即 $x=1$，$P_l(1)=1$ 对所有 $l$ 成立）处，取前 3 项得 $\frac{3}{2} - \frac{7}{8} + \frac{11}{16} = \frac{21}{16} = 1.3125$，而真值是 $1$。加到 $P_7$、$P_9$ 会继续往下压，但压得很慢——因为 $\mathrm{sign}(x)$ 在 $x=0$ 处跳变，系数衰减只有 $O(1/l)$ 量级。这就是后面要重点说的"系数衰减慢导致的端点收敛困难"。

5.4 多极展开其实是同一件事的另一种写法

把上面那套东西用在"远处看一个局域电荷分布"的场景，就得到多极展开。从 $\frac{1}{|\mathbf{r}-\mathbf{r}'|} = \frac{1}{r}\sum (r'/r)^l P_l(\cos\gamma)$ 出发，把电荷密度积进去，$l=0$ 项给单极矩（总电荷），$l=1$ 项给偶极矩，$l=2$ 项给四极矩。每一项的角向依赖恰好就是 $P_l(\cos\gamma)$。

这里有个实用心得：判断多极展开需要保留到第几阶，看的是源区尺寸与观察距离的比值 $r'/r$。这个比值每提高一阶就多乘一次，如果 $r'/r = 0.1$，那 $l=3$ 之后的项贡献已经在 $10^{-4}$ 以下，算到四极矩完全够用。反过来，如果是近场问题（$r' \approx r$），这个级数收敛极慢，老老实实做数值积分比展开更划算。

6. 最佳平方逼近：投影就是最优解

6.1 正规方程在正交基下的退化

回到拟合问题。给定 $f(x)$ 在 $[-1,1]$ 上，用 $\sum_{n=0}^{N} c_n P_n(x)$ 去逼近它，要求最小化平方误差：

$$E = \int_{-1}^{1}\left[f(x) - \sum_{n=0}^{N} c_n P_n(x)\right]^2 dx$$

对每个 $c_k$ 求偏导并令其为零，得到正规方程 $\sum_n c_n \langle P_n, P_k\rangle = \langle f, P_k\rangle$。由于 $\langle P_n, P_k\rangle = \frac{2}{2k+1}\delta_{nk}$，这个方程组直接塌缩成 $N+1$ 个标量等式：

$$c_k = \frac{2k+1}{2}\int_{-1}^{1} f(x)P_k(x),dx$$

这就是截断展开的系数就是最佳平方逼近系数的证明。换句话说，你不需要做任何优化迭代，直接把这个积分算出来，得到的就是 $L^2$ 意义下的最优解。这一点在数学上非常干净：正交投影就是最佳逼近，没有例外。

6.2 $|x|$ 与 $e^x$ 的收敛速度实测

我做过一组对比，结果很能说明问题。两个目标函数都是标准测试函数，但收敛速度差了三个数量级：

原因很直接：勒让德级数的收敛速度取决于函数的解析性。$e^x$ 在整个复平面上解析，系数按 $1/n!$ 衰减，十几项就能压到机器精度。$|x|$ 在 $x=0$ 处一阶导数不连续，收敛速度退化到幂律，需要几百项。$\mathrm{sign}(x)$ 直接有跳跃间断，在跳变点附近永远收敛不了，还会出现类似吉布斯现象的过冲。

实测数据：$|x|$ 的展开，$c_0 = 1/2$，$c_2 = 5/8$，之后偶数项按 $O(n^{-2})$ 衰减，$N=60$ 时 $L^2$ 误差约 $10^{-3}$，$N=600$ 才进入 $10^{-5}$。这个量级关系值得记住：只要目标函数存在一处导数不连续，展开项数就得按平方关系往上加。

6.3 与泰勒展开、切比雪夫展开的横向对比

勒让德不是唯一选择，实际工程里经常要在三套方案间挑：

经验判断：谱方法解微分方程优先切比雪夫（因为要控制最大误差），球坐标/轴对称物理问题优先勒让德（因为和几何结构天然对齐），需要符号运算或理论推导时看泰勒（因为形式最简洁）。有人试过在球坐标问题里硬用切比雪夫，结果边界条件匹配那一步就得重新推导一遍正交关系，得不偿失。

7. 高斯-勒让德求积：节点取零点的深层原因

7.1 $2n-1$ 次精度的构造性证明

$n$ 个节点最多能让求积公式精确到多少阶多项式？插值型求积的答案是 $n-1$，但高斯-勒让德把节点位置也当成自由度一起优化，直接干到 $2n-1$。

证明思路很短但很关键。取节点 $x_i$ 为 $P_n$ 的 $n$ 个零点，对任意次数不超过 $2n-1$ 的多项式 $f$，做带余除法：

$$f(x) = q(x)P_n(x) + r(x)$$

其中 $\deg q \le n-1$，$\deg r \le n-1$。两边积分，第一项 $\int_{-1}^{1} q(x)P_n(x)dx = 0$——因为 $q$ 的次数不超过 $n-1$，而 $P_n$ 与所有次数低于 $n$ 的多项式正交。所以 $\int f = \int r$。而 $r$ 的次数不超过 $n-1$，$n$ 点插值型求积对它精确。证毕。

这个论证的美妙之处在于：正交性直接换来了两倍的代数精度。同样是 $n$ 个函数求值，辛普森法则只有 4 阶精度，高斯-勒让德有 $2n-1$ 阶。

7.2 节点与权重的数值计算

权重公式是：

$$w_i = \frac{2}{(1-x_i^2)\left[P_n'(x_i)\right]^2}$$

手算小规模时常用这几组值：

注意权重之和恒等于 2，这正好是 $\int_{-1}^{1} 1,dx$。写代码时这是个现成的自检项：如果权重和不是 2，节点或权重的计算一定出了问题。

实际使用中，节点和权重一般用numpy.polynomial.legendre.leggauss(n)直接取。它的底层是用牛顿迭代求 $P_n$ 的零点，初值取 $\cos\left(\frac{(4i-1)\pi}{4n+2}\right)$，通常 3 到 5 次迭代就能收敛到机器精度。$n$ 到 100 左右都没有明显精度损失；$n$ 再往上，节点在端点附近会挤得很密，权重变得很小，求和时的相对误差开始显现。

7.3 区间映射和几个常见误用

高斯-勒让德的标准形式定义在 $[-1,1]$ 上。要算 $\int_a^b f(t)dt$，必须做映射 $t = \frac{b-a}{2}x + \frac{a+b}{2}$，并且乘上雅可比 $\frac{b-a}{2}$。漏掉这个因子是最常见的错误，表现形式是所有积分结果都差一个常数倍。

另外两个坑值得单独提：

• 节点位置被当成等距：很多人潜意识里觉得"节点应该是均匀的"，实际高斯节点在端点附近密集、中间稀疏。如果画出节点分布图发现是均匀的，说明拿到的是牛顿-科茨公式的节点，不是高斯节点。
• 被积函数有奇点却直接用：高斯-勒让德不处理端点奇点。$\int_0^1 x^{-1/2}dx$ 这种，直接跑会得到严重偏差。正确做法是做变量替换消掉奇点，比如取 $x = u^2$，$dx = 2u,du$，被积函数变成 $2$，光滑可积。

我实测过一个对比：$\int_0^1 e^{-x^2}dx$ 用 10 点高斯-勒让德，结果与参考值差 $10^{-16}$ 量级，几乎就是机器精度；同样 10 点的复合辛普森法则误差在 $10^{-10}$ 附近。点数越多，差距越明显。

8. 手写一份实现：递推核与三组自检

8.1 递推核的向量化写法

下面是能直接跑的完整实现，重点在于向量化处理和自检项：

import numpy as np def legendre_vals(n_max, x): """返回 P_0..P_{n_max} 在 x 处的值，形状 (n_max+1, len(x))""" x = np.asarray(x, dtype=float) P = np.empty((n_max + 1, x.size)) P[0] = 1.0 if n_max == 0: return P P[1] = x for k in range(1, n_max): P[k + 1] = ((2 * k + 1) * x * P[k] - k * P[k - 1]) / (k + 1) return P def legendre_coeffs(f, n_max, a=-1.0, b=1.0, n_quad=400): """计算 f 在 [a,b] 上的勒让德展开系数（含雅可比）""" xg, wg = np.polynomial.legendre.leggauss(n_quad) t = 0.5 * (b - a) * xg + 0.5 * (a + b) # 映射回物理区间 jac = 0.5 * (b - a) # 别忘了这个 ft = f(t) P = legendre_vals(n_max, xg) c = np.array([ (2 * k + 1) / 2.0 * jac * np.sum(wg * ft * P[k]) for k in range(n_max + 1) ]) return c

递推循环里每一步只依赖前两项，P数组预先分配好，避免在循环里追加，这样 $n=200$ 加 1000 个采样点也就毫秒级。系数计算里我把雅可比jac单独拎出来，就是为了让它在代码审查时一眼可见。

8.2 三组自检：正交性、重构、积分精度

写完核心函数后，我一定跑这三组检查，任何一组不过就回头查：

第一组，正交性检查。数值算 $\int_{-1}^{1}P_m P_n dx$ 与理论值 $\frac{2}{2n+1}\delta_{mn}$ 对比：

xg, wg = np.polynomial.legendre.leggauss(200) P = legendre_vals(15, xg) G = np.array([[np.sum(wg * P[i] * P[j]) for j in range(16)] for i in range(16)]) D = np.diag(2.0 / (2 * np.arange(16) + 1.0)) err = np.abs(G - D).max() # 预期 < 1e-13

跑下来误差通常在 $10^{-15}$ 量级。如果跑到 $10^{-8}$，说明递推写错了项系数，或者是浮点累积问题。

第二组，重构检查。取 $f(x) = e^x$，展开到 12 阶，在 1000 个点上重构，看最大偏差：

f = np.exp c = legendre_coeffs(f, 12) xt = np.linspace(-1, 1, 1000) approx = c @ legendre_vals(12, xt) print(np.abs(approx - f(xt)).max()) # 预期 ~1e-15

第三组，积分精度检查。用 8 点高斯-勒让德算 $\int_{-1}^{1}x^{14}dx$，理论值 $\frac{2}{15}$。8 点可精确到 15 次，14 次在射程内，误差应在 $10^{-15}$。

这三组检查覆盖了正交性、展开、积分三条主线，加起来不到 30 行代码，每次改完实现都值得重跑一遍。

8.3 高阶场合的稳定性处理

阶数上到 100 以上时，有两件事要注意。一是函数值的量级：虽然理论上 $|P_n|\le 1$，但中间递推过程里的乘法和累加会带来舍入误差积累，实测 $n=300$ 时误差在 $10^{-11}$ 量级，仍然可以接受。二是系数的动态范围：如果函数在端点附近变化剧烈，$c_n$ 可能从 $10^{0}$ 一路降到 $10^{-16}$，此时直接求和会因为大数吃小数丢失精度。稳妥的做法是按阶数分块求和，或者用math.fsum做精确累加。

还有一个更隐蔽的问题：递推方向。上面用的是从低阶到高阶的递推，这在 $|x| \le 1$ 时是稳定的。但如果有人图省事从高阶往下推（用同一个式子反解 $P_{n-1}$），误差会指数放大，$n=50$ 就完全不可用。这个方向性来源于递推关系的两个特征解一个递增一个递减，数值上必须顺着增长的方向推。

9. 几个真踩过的坑和对应的处理方式

9.1 "系数很小"不代表"误差很小"

这是我栽得最深的一次。做某组球面数据的重构时，看到 $c_{20}$ 已经降到 $10^{-6}$ 就判定收敛了，直接截断在 20 阶。结果在 $\theta$ 接近极点的地方，重构值和原数据差到了 $10^{-2}$。回头一查才明白：系数的衰减速度只控制整体的 $L^2$ 误差，不控制逐点误差。

原因在上一节那个双半球算例里已经露头了。端点上 $P_l(1)=1$ 对所有 $l$ 成立，意味着每一项在端点都是满额贡献，逐项相加时是交替级数，收敛速度取决于系数衰减，而系数衰减只有 $O(1/l)$。$l=20$ 时截断误差理论上还在 $1/20$ 量级附近，虽然 $c_{20}$ 本身很小，但后面的项加起来仍然可观。

正确的做法是：除了看系数，必须做逐点误差检查。我现在的习惯是在原始数据密集的地方取 200 个点，算 $|f - \text{重构}|$ 的最大值，这个数达标了才算收敛。只看系数会被骗得很惨。

9.2 区间映射的雅可比漏乘，图形上看不出来

前面提过一次，这里展开说。漏乘雅可比的典型症状是"形状完美、幅值不对"。因为所有系数都差同一个比例因子，重构曲线和目标函数在归一化之后完全重合，眼看看不出问题。

定位方法很简单：在重构结果上除一个常数做缩放，看能不能对上。如果缩放后完全重合，那就是漏了雅可比或者归一化因子；如果缩放后仍对不上，说明是基函数的构造有问题。我一般直接打印前三个系数，和手算值对比，手算三个数比调试半天快得多。

9.3 采样点数与求积阶数的隐性耦合

计算系数时我用高斯-勒让德积分去算那个 $\int f P_k,dx$。这里有个容易忽略的点：求积节点数必须显著大于展开阶数。如果展开到 20 阶，求积只用 20 个点，那对 $f \cdot P_{20}$ 这个 20 阶以上多项式，20 点高斯积分的精度就未必够。

经验规则是求积节点数取展开阶数的 2 到 4 倍，稳妥一些。代价不大，但能避免一类很难查的误差——表现为高阶系数出现不规则的抖动，而不是平滑衰减。

9.4 高阶展开时的条件数陷阱

有个反直觉的现象：正交基解决了正规方程的病态，但如果直接用 $A^{\mathsf T}A$ 求逆的方式算系数（而不是用正交性逐个积分），高阶时仍然会出问题，因为 $A^{\mathsf T}A$ 虽然在理论上是对角阵，数值上仍有 $10^{-14}$ 量级的非对角残余，求逆会把它放大。

所以永远走积分公式那条路，别走矩阵求逆。即便形式上都等价，数值行为完全不同。这是我在两套实现上对比出来的：同样 30 阶，积分法误差 $10^{-14}$，矩阵求逆误差 $10^{-5}$，差了九个数量级。

9.5 关于"该用几阶"的经验判断

最后分享一个判断流程，能省掉大量试错。先看目标函数在区间内有没有导数不连续点：有的话，先绕开它做分段展开，一般比硬堆阶数划算得多。没有的话，看特征尺度，函数变化一个完整周期大概需要 8 到 12 项，拿"区间宽度除以最小特征长度"再乘 10 就是起始阶数。然后逐步加阶，同时监控逐点误差和系数衰减速度，系数衰减到 $10^{-13}$ 左右基本就是双精度的天花板了，再往上加没有意义。

对于物理问题，还有个更快的捷径：从问题的几何结构直接判断。如果是一个平滑场，主要贡献集中在低阶多极矩，$l=10$ 以内通常够；如果有尖锐边界或局部源，就得准备 $l=50$ 往上的规模，并且提前为端点附近的收敛问题做好心理准备。