基于非等距切比雪夫节点的重心权公式数值微分及其稳定性分析
题目描述
在数值微分中，我们通常需要根据给定的一组离散函数值 \(\{ (x_i, f(x_i)) \}_{i=0}^{n}\) 来近似计算函数在某个点（可能是节点，也可能是内点）处的导数值。对于高精度要求，特别是当需要在整个区间上获得导数近似时，我们常采用基于多项式插值的微分公式。当节点为高精度的插值节点（如切比雪夫节点）时，可以得到更稳定的近似。重心权公式是计算拉格朗日插值及其导数的一种高效且数值稳定的方法。本题将引导你理解并构造基于非等距节点（具体为切比雪夫节点）的数值微分公式，重点阐述重心权公式的构建方法，并分析其相对于标准拉格朗日插值微分公式的稳定性优势。
解题过程
步骤1：问题设定与切比雪夫节点选取
我们的目标是在区间 \([-1, 1]\) 上近似计算一个光滑函数 \(f(x)\) 的导数 \(f'(x)\)。我们并非在等距节点上采样，而是采用第一类切比雪夫节点，因为它们能最小化龙格（Runge）现象，并通常能带来更高的插值精度和更好的数值稳定性。第 \(k\) 个 \(n+1\) 个节点的切比雪夫节点定义为：
\[x_k = \cos\left(\frac{2k+1}{2n+2} \pi\right), \quad k = 0, 1, \dots, n \]
这些节点是区间 \([-1,1]\) 内的点，在端点处更密集。
给定函数在这些节点上的值 \(f_k = f(x_k)\)，我们希望计算在任意点 \(z \in [-1,1]\) 处（尤其可能在节点上）的导数近似值 \(f'(z) \approx p_n'(z)\)，其中 \(p_n(x)\) 是经过这 \(n+1\) 个节点的 \(n\) 次拉格朗日插值多项式。
步骤2：标准拉格朗日插值微分公式及其潜在问题
拉格朗日插值多项式为：
\[p_n(x) = \sum_{j=0}^{n} f_j L_j(x) \]
其中 \(L_j(x) = \prod_{\substack{i=0 \\ i \neq j}}^{n} \frac{x - x_i}{x_j - x_i}\) 是拉格朗日基多项式。
在点 \(z\) 处的导数近似为：
\[p_n'(z) = \sum_{j=0}^{n} f_j L_j'(z) \]
若 \(z = x_k\) 是一个节点，则有节点微分公式：
\[f'(x_k) \approx p_n'(x_k) = \sum_{j=0}^{n} f_j L_j'(x_k) = \sum_{j=0}^{n} D_{kj} f_j \]
其中 \(D_{kj} = L_j'(x_k)\) 构成了微分矩阵。
直接计算基函数的导数公式为：
\[L_j'(z) = L_j(z) \sum_{\substack{i=0 \\ i \neq j}}^{n} \frac{1}{z - x_i} \]
而当 \(z = x_k\) 时，有
\[D_{kj} = L_j'(x_k) = \begin{cases} \frac{1}{x_k - x_j} \prod_{\substack{i=0 \\ i \neq j, k}}^{n} \frac{x_k - x_i}{x_j - x_i}, & j \neq k \\ \sum_{\substack{i=0 \\ i \neq k}}^{n} \frac{1}{x_k - x_i}, & j = k \end{cases} \]
潜在问题：当节点数 \(n\) 较大时，直接计算上述乘积和商可能会由于舍入误差的积累而变得不稳定，尤其对于等距节点。即使对于切比雪夫节点，虽然插值更稳定，但直接使用上述显式公式计算 \(D_{kj}\) 仍可能引入不必要的舍入误差。
步骤3：重心权公式的引入
重心权公式为拉格朗日插值提供了一个在数值上更稳定的表达方式。首先定义重心权重 \(w_j\)：
\[w_j = \frac{1}{\prod_{\substack{i=0 \\ i \neq j}}^{n} (x_j - x_i)}, \quad j = 0, \dots, n \]
对于切比雪夫节点，其重心权重有一个简洁的解析形式：
\[w_j = (-1)^j \sin\left(\frac{2j+1}{2n+2} \pi\right) \]
（准确地说，常包含一个归一化因子，但通常我们使用上述带符号的权重，因为比例常数在比值中会消去）。
利用这些权重，插值多项式可以写成重心形式：
\[p_n(x) = \frac{\sum_{j=0}^{n} \frac{w_j}{x - x_j} f_j}{\sum_{j=0}^{n} \frac{w_j}{x - x_j}} \]
这个形式的优势在于，权重 \(w_j\) 只依赖于节点，不依赖于函数值，可以一次性预先计算并存储，且对切比雪夫节点是数值稳定的。
步骤4：基于重心形式的数值微分公式推导
我们的目标是计算 \(p_n'(z)\)。设：
\[u(x) = \sum_{j=0}^{n} \frac{w_j}{x - x_j} f_j, \quad v(x) = \sum_{j=0}^{n} \frac{w_j}{x - x_j} \]
则 \(p_n(x) = u(x) / v(x)\)。
利用商的求导法则：
\[p_n'(z) = \frac{u'(z) v(z) - u(z) v'(z)}{[v(z)]^2} \]
我们需要计算 \(u'(z)\), \(v'(z)\), \(u(z)\), \(v(z)\)。
计算 \(u(z)\) 和 \(v(z)\) 是直接的：
\[u(z) = \sum_{j=0}^{n} \frac{w_j}{z - x_j} f_j, \quad v(z) = \sum_{j=0}^{n} \frac{w_j}{z - x_j} \]
对于导数项，注意到：
\[\frac{d}{dx} \left( \frac{w_j}{x - x_j} \right) = -\frac{w_j}{(x - x_j)^2} \]
因此：
\[u'(z) = -\sum_{j=0}^{n} \frac{w_j}{(z - x_j)^2} f_j, \quad v'(z) = -\sum_{j=0}^{n} \frac{w_j}{(z - x_j)^2} \]
步骤5：简化与最终的通用公式
将以上结果代入求导公式，得到在任意点 \(z\) 处的导数近似：
\[p_n'(z) = \frac{ \left[ -\sum_{j=0}^{n} \frac{w_j}{(z - x_j)^2} f_j \right] \left[ \sum_{j=0}^{n} \frac{w_j}{z - x_j} \right] - \left[ \sum_{j=0}^{n} \frac{w_j}{z - x_j} f_j \right] \left[ -\sum_{j=0}^{n} \frac{w_j}{(z - x_j)^2} \right] }{ \left[ \sum_{j=0}^{n} \frac{w_j}{z - x_j} \right]^2 } \]
这可以进一步化简整理为：
\[p_n'(z) = \frac{ \sum_{j=0}^{n} \frac{w_j}{z - x_j} \left( f_j - p_n(z) \right) \cdot \left( \sum_{\substack{k=0 \\ k \neq j}}^{n} \frac{w_k / (z - x_k)}{w_j / (z - x_j)} \cdot \frac{1}{x_j - x_k}? \right) }{ \sum_{j=0}^{n} \frac{w_j}{z - x_j} } \quad \text{(此形式复杂，不常用)} \]
实际上，一个更实用且数值稳定的计算方式是直接使用商的求导法则计算出来的原始形式，但可以避免直接计算两个大数相减。一个常用的稳定实现方法是：
- 预先计算并存储重心权重 \(w_j\)。
- 对于给定的 \(z\)，计算：
\[ s_1 = \sum_{j=0}^{n} \frac{w_j}{z - x_j}, \quad s_2 = \sum_{j=0}^{n} \frac{w_j}{(z - x_j)^2} \]
- 计算插值函数值在 \(z\) 处的近似：\(p_z = \left( \sum_{j=0}^{n} \frac{w_j}{z - x_j} f_j \right) / s_1\)。
- 然后导数近似为：
\[ p_n'(z) = \frac{ \left( \sum_{j=0}^{n} \frac{w_j}{z - x_j} f_j \right) \cdot s_2 / s_1 - \sum_{j=0}^{n} \frac{w_j}{(z - x_j)^2} f_j }{ s_1 } \]
为了避免计算 \(p_z\) 的误差传播，可以重新排列：
\[ p_n'(z) = \frac{1}{s_1} \left( \frac{ \left( \sum_{j=0}^{n} \frac{w_j}{z - x_j} f_j \right) s_2 }{ s_1 } - \sum_{j=0}^{n} \frac{w_j}{(z - x_j)^2} f_j \right) \]
但更常见且稳定的方式是直接计算：
\[ p_n'(z) = \frac{ \sum_{j=0}^{n} \frac{w_j}{z - x_j} (f_j - p_z) \cdot \frac{1}{z - x_j}? }{ s_1 } \quad \text{(不精确)} \]
实际上，标准且稳定的公式是：
\[ p_n'(z) = \frac{ \sum_{j=0}^{n} \frac{w_j}{z - x_j} \left( \frac{f_j - p_z}{z - x_j} \right) }{ s_1 } \]
验证：从 \(p_n(z) = p_z\) 出发，考虑基函数表达式的导数公式的稳定形式。一个已知的稳定公式是：
\[ p_n'(z) = \frac{ \sum_{j=0}^{n} \frac{w_j}{z - x_j} \left( \frac{f_j - p_n(z)}{z - x_j} \right) }{ \sum_{j=0}^{n} \frac{w_j}{z - x_j} } \]
这个形式避免了直接计算 \(u'(z)v(z) - u(z)v'(z)\) 时可能出现的相近大数相减，是推荐的计算方式。
步骤6：节点处的导数公式简化
当 \(z = x_k\) 时，上述公式存在 \(0/0\) 型未定式，需要取极限。此时，稳定公式简化为：
\[p_n'(x_k) = \frac{ \frac{w_k'}{w_k} f_k + \sum_{\substack{j=0 \\ j \neq k}}^{n} \frac{w_j}{x_k - x_j} \left( \frac{f_j - f_k}{x_k - x_j} \right) }{ \frac{w_k'}{w_k} + \sum_{\substack{j=0 \\ j \neq k}}^{n} \frac{w_j}{x_k - x_j} \cdot \frac{1}{x_k - x_j}? } \]
实际上，更简单的做法是直接使用微分矩阵 \(D_{kj}\)，但其元素可以用重心权重更稳定地计算：
对于 \(j \neq k\)：
\[D_{kj} = \frac{w_j / w_k}{x_k - x_j} \]
对于对角线元素 \(j = k\)：
\[D_{kk} = -\sum_{\substack{j=0 \\ j \neq k}}^{n} D_{kj} \]
这个计算 \(D_{kj}\) 的方法比直接乘积公式在数值上更稳定，因为 \(w_j\) 的计算（特别是对于切比雪夫节点）是良态的。
步骤7：稳定性分析
-
重心权公式的优势：
- 数值稳定：标准拉格朗日插值在计算基函数或其导数时，涉及大量乘除法，容易导致上溢/下溢和舍入误差累积。重心形式将大部分计算转化为对预计算的权重 \(w_j\) 的加权求和，减少了运算次数和舍入误差。
- 高效：一旦节点确定，权重 \(w_j\) 只需计算一次并存储。对于不同的 \(z\) 或不同的函数 \(f\)，都可以重用这些权重。
- 适用于切比雪夫节点：切比雪夫节点的重心权重有简洁的解析表达式，计算快速准确，且权重的数量级变化相对温和，进一步增强了稳定性。
-
与等距节点的对比：等距节点的重心权重随着 \(n\) 增大变化剧烈（例如，在端点附近权重绝对值非常大），容易导致严重的舍入误差。而切比雪夫节点的权重变化平缓，使得基于它们的重心权公式在整个区间上都能保持高精度。
-
计算 \(p_n'(z)\) 的稳定性：推荐的公式 \(p_n'(z) = \frac{ \sum_{j=0}^{n} \frac{w_j}{z - x_j} \left( \frac{f_j - p_n(z)}{z - x_j} \right) }{ \sum_{j=0}^{n} \frac{w_j}{z - x_j} }\) 通过先计算插值值 \(p_n(z)\)，然后在求和时减去它，避免了分子中两个大数（\(u'(z)v(z)\) 和 \(u(z)v'(z)\)）的直接相减，这是数值计算中提高精度的常用技巧。
-
适用性与限制：该方法适用于任意分布节点，但对切比雪夫节点特别有效。当 \(n\) 非常大时（例如几百），任何基于高次多项式插值的微分都会对函数值的微小误差（噪声）极度敏感，这是多项式插值的固有特性，并非重心形式特有。但对于光滑函数和适度大小的 \(n\)，基于切比雪夫节点和重心权公式的数值微分是稳定且高精度的。
总结：本题展示了如何利用非等距的切比雪夫节点和重心权公式来构造稳定、高效的数值微分方法。核心在于通过重心权重重新组织计算流程，避免了直接处理拉格朗日基函数的高阶多项式形式，从而在数值实现上显著提高了稳定性，尤其适用于需要在整个区间上获取导数近似的高精度计算场景。