基于切比雪夫展开的快速振荡积分计算：Clenshaw-Curtis 积分法
题目描述
计算定积分
\[I = \int_{-1}^{1} f(x) \, dx \]
其中，被积函数 $ f(x) $ 是一个光滑函数，但在区间 $[-1, 1]$ 上可能表现出快速振荡特性，或者其高阶导数变化剧烈，使得传统的牛顿-柯特斯公式（如辛普森法）需要极细的剖分才能达到所需精度，计算效率低下。要求设计一种高效、高精度的数值积分方法。
解题过程循序渐进讲解
我们解决此问题的核心思路是：将函数用一组特别适合数值计算的正交多项式（切比雪夫多项式）展开，然后利用展开式系数的递推关系，高效、高精度地计算积分的近似值。 这种方法称为 Clenshaw-Curtis 积分法。
第一步：理解理论基础——切比雪夫多项式与函数逼近
-
切比雪夫多项式：定义在区间 $[-1, 1]$ 上的一组重要多项式，通常记为 $ T_n(x) $，其中 $ n $ 是多项式的次数。
- 其定义为：$ T_n(x) = \cos(n \arccos x) $。这个定义巧妙地将多项式与三角函数联系起来。
- 前几项为：$ T_0(x) = 1 $， $ T_1(x) = x $， $ T_2(x) = 2x^2 - 1 $，等等。
- 重要性质：它们是关于权函数 $ w(x) = 1/\sqrt{1 - x^2} $ 正交的多项式。但在 Clenshaw-Curtis 积分中，我们更关键地利用它们的离散正交性和与傅里叶级数的紧密联系。
-
函数展开：任何在 $[-1, 1]$ 上平方可积的函数 $ f(x) $ 都可以展开成切比雪夫级数：
\[ f(x) = \frac{a_0}{2} + \sum_{k=1}^{\infty} a_k T_k(x) \]
理论上，展开系数 \$ a_k \$ 由内积 \$ \langle f, T_k \rangle \$ 给出。但在数值计算中，我们无法计算无穷级数，需要进行截断。
第二步：从连续展开到离散采样——切比雪夫插值
为了计算系数 $ a_k $，Clenshaw-Curtis 积分法采用了一个关键策略：在切比雪夫点上对函数进行采样，然后通过离散余弦变换 (DCT) 来近似计算展开系数。
- 切比雪夫点：在区间 $[-1, 1]$ 上，$ N+1 $ 个切比雪夫点定义为：
\[ x_j = \cos\left( \frac{j\pi}{N} \right), \quad j = 0, 1, \dots, N \]
这些点不是均匀分布的，而是在区间端点处更密集。这种分布能有效抑制高次多项式插值在区间端点附近的龙格现象，使得插值过程数值稳定。
- 离散插值与系数近似：我们构造一个 $ N $ 次多项式 $ p_N(x) $，使其在切比雪夫点 $ x_j $ 处精确等于 $ f(x_j) $。可以证明，这个插值多项式可以写为：
\[ p_N(x) = \sum_{k=0}^{N} {}' c_k T_k(x) \]
这里的 \$ \sum{}' \$ 表示求和时第一项（\$ k=0 \$）的系数要除以 2。系数 \$ c_k \$ 是**离散切比雪夫系数**，它是对理论系数 \$ a_k \$ 的近似。最关键的是，\$ c_k \$ 可以通过函数在切比雪夫点处的采样值 \$ f(x_j) \$ 高效计算：
\[ c_k = \frac{2}{N} \sum_{j=0}^{N} {}'' f(x_j) \cos\left( \frac{k j \pi}{N} \right) \]
这里的 \$ \sum{}'' \$ 表示求和时第一项（\$ j=0 \$）和最后一项（\$ j=N \$）的权重为 1/2。这个计算形式正是**离散余弦变换 (DCT-I)**，存在快速算法 (FFT)，计算复杂度为 \$ O(N \log N) \$，非常高效。
第三步：从展开式到积分近似
现在，我们用插值多项式 $ p_N(x) $ 来近似原函数 $ f(x) $，并用 $ p_N(x) $ 的积分来近似原积分 $ I $。
\[I = \int_{-1}^{1} f(x) dx \approx \int_{-1}^{1} p_N(x) dx = \int_{-1}^{1} \left[ \sum_{k=0}^{N} {}' c_k T_k(x) \right] dx \]
利用积分的线性性质，交换积分与求和顺序：
\[I \approx \sum_{k=0}^{N} {}' c_k \int_{-1}^{1} T_k(x) dx \]
第四步：计算切比雪夫多项式的积分——递推公式
我们需要计算积分 $ \int_{-1}^{1} T_k(x) dx $。利用切比雪夫多项式的定义和三角恒等式，可以得到一个简洁的结果：
- 当 $ k = 0 $ 时，$ \int_{-1}^{1} T_0(x) dx = \int_{-1}^{1} 1 dx = 2$。
- 当 $ k = 1 $ 时，$ \int_{-1}^{1} T_1(x) dx = \int_{-1}^{1} x dx = 0$。
- 对于 $ k \ge 2 $，利用递推关系 $ T_{k+1}(x) = 2x T_k(x) - T_{k-1}(x) $ 和分部积分，可以推导出：
\[ \int_{-1}^{1} T_k(x) dx = \begin{cases} \frac{2}{1-k^2}, & \text{$k$ 为偶数} \\ 0, & \text{$k$ 为奇数} \end{cases} \]
注意，分母是 \$ 1 - k^2 \$，当 \$ k=1 \$ 时也为0，与上述结果一致。
第五步：得到 Clenshaw-Curtis 积分公式
将积分结果代入第三步的近似式：
\[I \approx \sum_{k=0}^{N} {}' c_k \cdot w_k \]
其中，权重 $ w_k = \int_{-1}^{1} T_k(x) dx $，即：
\[w_k = \begin{cases} 2, & k=0 \\ 0, & k \text{ 为奇数} \\ \frac{2}{1-k^2}, & k \text{ 为偶数}, k \ge 2 \end{cases} \]
由于当 $ k $ 为奇数时权重为零，求和实际上只在偶数 $ k $ 上进行。最终公式为：
\[I \approx c_0 \cdot 1 + \sum_{\substack{k=2 \\ k \text{ even}}}^{N} c_k \cdot \frac{2}{1-k^2} \]
注意，这里利用了 $ \sum{}' $ 的约定，$ c_0 $ 项实际上对应于 $ (c_0/2) * 2 = c_0 $。
第六步：算法步骤总结与特性
- 选择点数 N：通常选择 $ N = 2^m $（如 16, 32, 64, ...），以便利用 FFT 加速 DCT 计算。
- 生成采样点：计算切比雪夫点 $ x_j = \cos(j\pi/N), j=0,...,N$。
- 函数求值：计算函数在这些点上的值 $ f_j = f(x_j)$。
- 计算离散切比雪夫系数：通过 DCT-I 计算系数 $ c_k$：
\[ c_k = \frac{2}{N} \left[ \frac{f_0}{2} + \sum_{j=1}^{N-1} f_j \cos\left( \frac{k j \pi}{N} \right) + \frac{f_N}{2} \cos(k\pi) \right] \]
- 计算积分近似值：
\[ I_N = c_0 + \sum_{\substack{k=2 \\ k \text{ even}}}^{N} \frac{2c_k}{1-k^2} \]
- 精度控制（可选自适应）：可以比较不同 $ N $（例如 $ N $ 和 $ 2N $）下的结果 $ I_N $ 和 $ I_{2N} $。如果其差值的绝对值小于预设容差，则接受 $ I_{2N} $ 为最终结果；否则，将 $ N $ 加倍，重复步骤 2-5。
方法特性：
- 高精度与指数收敛：如果 $ f(x) $ 是解析函数或非常光滑，Clenshaw-Curtis 积分的误差通常以 $ O(\rho^{-N}) $ 的速度指数衰减（$ \rho > 1 $），收敛速度远快于多项式收敛的牛顿-柯特斯公式。
- 高效稳定：核心计算是 DCT，可用 FFT 在 $ O(N \log N) $ 时间内完成，且切比雪夫点采样避免了高次多项式插值的不稳定性。
- 易于嵌套：当点数 $ N $ 加倍时，新的切比雪夫点集合完全包含旧的点集，这意味着所有旧的函数采样值都可以复用，只需计算新增点的函数值。这一特性使其非常适合实现自适应积分，在需要更高精度的地方（通过加密点数）能高效利用已有计算。
通过以上步骤，我们就将一个困难的振荡函数积分问题，转化为了在特殊点集上采样、进行快速变换、然后加权求和的稳定高效过程，这就是 Clenshaw-Curtis 积分法的精髓。