基于 Clenshaw–Curtis 积分法的高效节点与权重预计算及其在傅里叶系数递推中的应用
题目描述
考虑数值计算积分
\[I = \int_{-1}^{1} f(x) \, dx. \]
Clenshaw–Curtis 求积公式是一种基于切比雪夫节点（即 \(x_k = \cos(k\pi/N)\)，\(k=0,1,\dots,N\)）的插值型求积公式。本题目聚焦于该公式的两个核心高效实现技巧：
- 如何利用快速傅里叶变换（FFT）高效计算求积权重，避免直接求解线性方程组。
- 如何利用切比雪夫多项式的递推关系，将函数值的加权和表达为傅里叶系数的线性组合，从而通过一次 FFT 同时得到所有节点处的函数值和积分近似值。
请循序渐进地讲解其原理、推导过程与算法步骤。
解题过程
第一步：Clenshaw–Curtis 求积公式的基本形式
Clenshaw–Curtis 求积公式采用 \(N+1\) 个切比雪夫节点（第一类）：
\[x_k = \cos\left( \frac{k\pi}{N} \right), \quad k=0,1,\dots,N. \]
在这些节点上对 \(f(x)\) 进行拉格朗日插值，然后对插值多项式精确积分，得到求积公式：
\[I_N = \sum_{k=0}^{N} w_k f(x_k), \]
其中 \(w_k\) 是求积权重。直接通过拉格朗日基函数积分计算 \(w_k\) 的复杂度为 \(O(N^2)\)。下面介绍一种基于 FFT 的 \(O(N \log N)\) 高效算法。
第二步：利用切比雪夫展开表示被积函数
将 \(f(x)\) 用切比雪夫多项式 \(T_n(x) = \cos(n \arccos x)\) 展开：
\[f(x) \approx \sum_{n=0}^{N} {}' a_n T_n(x), \]
其中求和符号上的撇号表示首项系数减半，即：
\[\sum_{n=0}^{N} {}' a_n T_n(x) = \frac{a_0}{2} + \sum_{n=1}^{N-1} a_n T_n(x) + \frac{a_N}{2} T_N(x). \]
展开系数 \(a_n\) 由离散余弦变换（DCT）给出：
\[a_n = \frac{2}{N} \sum_{k=0}^{N} {}'' f(x_k) T_n(x_k) = \frac{2}{N} \sum_{k=0}^{N} {}'' f\left( \cos\frac{k\pi}{N} \right) \cos\frac{nk\pi}{N}, \]
这里双撇号表示求和时首末项减半，即：
\[\sum_{k=0}^{N} {}^{''} g_k = \frac{g_0}{2} + \sum_{k=1}^{N-1} g_k + \frac{g_N}{2}. \]
第三步：积分值的表达式
对切比雪夫展开式在 \([-1,1]\) 上逐项积分。利用切比雪夫多项式的正交性：
\[\int_{-1}^{1} T_n(x) \, dx = \begin{cases} 2, & n=0, \\ 0, & n \text{ 为奇数}, \\ \frac{2}{1-n^2}, & n \text{ 为偶数且 } n \ge 2. \end{cases} \]
因此，
\[I = \int_{-1}^{1} f(x) dx \approx \int_{-1}^{1} \sum_{n=0}^{N} {}' a_n T_n(x) \, dx = \sum_{n=0}^{N} {}' a_n \int_{-1}^{1} T_n(x) \, dx. \]
将积分值记为 \(b_n = \int_{-1}^{1} T_n(x) dx\)，则
\[I_N = \sum_{n=0}^{N} {}' a_n b_n. \]
代入 \(b_n\) 的表达式：
\[I_N = 2a_0 + \sum_{\substack{n=2 \\ n \text{ even}}}^{N} a_n \cdot \frac{2}{1-n^2}. \]
注意这里求和时 \(n\) 为偶数，且当 \(N\) 为偶数时，最后一项 \(n=N\) 的系数为 \(\frac{2}{1-N^2}\)（因为 \(a_N\) 已减半，但公式中已考虑“\('\)”符号，因此直接代入即可）。
第四步：通过 FFT 同时计算 \(a_n\) 和 \(I_N\)
- 计算函数值向量：在节点 \(x_k = \cos(k\pi/N)\) 上计算 \(f_k = f(x_k)\)，\(k=0,\dots,N\)。
- 应用 DCT：系数 \(a_n\) 可通过 DCT 得到。具体而言，定义向量 \(\mathbf{f} = (f_0, f_1, \dots, f_N)^T\)，则
\[ a_n = \frac{2}{N} \cdot \text{DCT}(\mathbf{f})_n, \quad n=0,\dots,N, \]
其中 DCT 可以通过 FFT 在 \(O(N \log N)\) 时间内计算。
3. 计算积分近似：
\[ I_N = 2a_0 + \sum_{\substack{m=1 \\ \text{或等价地处理偶数指标}}}^{\lfloor N/2 \rfloor} a_{2m} \cdot \frac{2}{1-(2m)^2}. \]
注意：实际计算时，通常预计算权重向量 \(\mathbf{v} = (v_0, v_1, \dots, v_N)\)，其中
\[ v_n = \begin{cases} 2, & n=0, \\ 0, & n \text{ 奇数}, \\ \frac{2}{1-n^2}, & n \text{ 偶数且 } n\ge 2, \end{cases} \]
则
\[ I_N = \sum_{n=0}^{N} {}' a_n v_n = \frac{a_0 v_0}{2} + \sum_{n=1}^{N-1} a_n v_n + \frac{a_N v_N}{2}. \]
由于 \(v_N\) 在 \(N\) 偶数时为 \(\frac{2}{1-N^2}\)，奇数时为 0，公式自动满足。
第五步：求积权重的显式表达（可选预计算）
虽然我们不需要显式计算权重即可得到积分值，但有时权重 \(w_k\) 本身也有用。可以证明：
\[w_k = \frac{2}{N} \sum_{n=0}^{N} {}'' \frac{1}{1-n^2} T_n(x_k) = \frac{2}{N} \left[ \frac{1}{2} + \sum_{n=1}^{N-1} \frac{\cos(n k\pi/N)}{1-n^2} + \frac{\cos(N k\pi/N)}{2(1-N^2)} \right], \]
其中最后一项仅当 \(N\) 为偶数时需要包含（因为当 \(N\) 为奇数时，\(1-N^2=0\) 导致该项为零权重分配，实际上公式中常将其归入求和方式处理）。
权重 \(w_k\) 也可以通过逆 DCT 快速计算，但通常我们只需用第四步的系数法计算积分，无需显式生成权重。
第六步：算法步骤总结
- 输入积分区间 \([-1,1]\)（如果是一般区间 \([a,b]\)，先做线性变换到 \([-1,1]\)），选择点数 \(N\)（通常为 2 的幂次加一，如 \(N=2^m+1\)）。
- 生成切比雪夫节点：
\[ x_k = \cos\left( \frac{k\pi}{N} \right), \quad k=0,\dots,N. \]
- 计算函数值 \(f_k = f(x_k)\)。
- 对序列 \(f_k\) 执行 DCT（通过 FFT 实现），得到系数 \(a_n\)：
\[ a_n = \frac{2}{N} \sum_{k=0}^{N} {}'' f_k \cos\left( \frac{nk\pi}{N} \right), \quad n=0,\dots,N. \]
- 预计算向量 \(v_n\)：
\[ v_0=2, \quad v_n=0 \ (\text{n 奇数}), \quad v_n = \frac{2}{1-n^2} \ (\text{n 偶数}, n\ge2). \]
- 计算积分近似：
\[ I_N = \sum_{n=0}^{N} {}' a_n v_n = \frac{a_0 v_0}{2} + \sum_{n=1}^{N-1} a_n v_n + \frac{a_N v_N}{2}. \]
当 \(N\) 为奇数时，\(v_N=0\)；当 \(N\) 为偶数时，\(v_N = \frac{2}{1-N^2}\)。
第七步：优点与注意事项
- 复杂度：主要开销在 DCT，可通过 FFT 在 \(O(N \log N)\) 时间完成，比直接高斯求积的节点生成（需解非线性方程）更高效，尤其当 \(N\) 较大且需重复计算多个相似积分时。
- 精度：Clenshaw–Curtis 求积具有多项式精度至少 \(N\) 阶，且当 \(f\) 光滑时误差以 \(O(N^{-k})\) 衰减，与高斯求积相近。
- 稳定性：权重均为正，数值稳定。
- 自适应扩展：可通过加倍节点数（如 \(N\) 取 2 的幂次加一）并复用已有函数值进行自适应加密，适合递归自适应实现。
通过以上步骤，我们不仅得到了 Clenshaw–Curtis 积分的快速算法，还理解了其背后基于切比雪夫展开和傅里叶技巧的数学原理。这种方法将数值积分转化为函数值的加权和，并通过快速变换极大提升了计算效率。