基于重心权公式的数值微分及其在非均匀节点下的高精度构造
题目描述
在科学计算中，函数导数经常无法解析求出，需用数值微分近似。中心差分公式在等距节点下简单有效，但当节点非均匀（如来自实验数据或自适应网格）时，传统差分公式精度下降或不适用。本题要求：推导基于非均匀节点拉格朗日插值的数值微分公式，特别是利用重心权公式（Barycentric Weights）高效构造高精度数值微分公式，分析其截断误差，并说明在非均匀节点下如何实现高精度微分逼近。
解题过程
第一步：问题分析与基本思路
数值微分的目标：给定函数 \(f(x)\) 在 \(n+1\) 个节点 \(x_0, x_1, \dots, x_n\)（未必等距）处的函数值 \(f_i = f(x_i)\)，求 \(f'(x)\) 在某个点（如节点或中间点）的近似值。
基本思路：先构造 \(f(x)\) 的插值多项式 \(p(x)\)，然后对 \(p(x)\) 求导作为 \(f'(x)\) 的近似。拉格朗日插值形式便于理论分析，但直接对其求导计算量较大（\(O(n^2)\) 次运算）。重心权公式是拉格朗日插值的一种等价但更稳定的表示，可高效计算插值及导数。
第二步：拉格朗日插值多项式及其导数形式
给定节点 \(\{x_i\}_{i=0}^n\)（互异），拉格朗日基函数：
\[\ell_i(x) = \prod_{\substack{j=0 \\ j \neq i}}^n \frac{x - x_j}{x_i - x_j}, \quad i=0,\dots,n. \]
插值多项式：
\[p(x) = \sum_{i=0}^n f_i \, \ell_i(x). \]
在节点 \(x_k\) 处对 \(p(x)\) 求导：
\[p'(x_k) = \sum_{i=0}^n f_i \, \ell_i'(x_k). \]
记 \(D_{ki} = \ell_i'(x_k)\)，则数值微分公式为线性组合：
\[f'(x_k) \approx \sum_{i=0}^n D_{ki} f_i. \]
关键是如何高效、稳定地计算 \(D_{ki}\)。
第三步：引入重心权与重心插值公式
定义重心权：
\[w_i = \frac{1}{\prod_{j \neq i} (x_i - x_j)}, \quad i=0,\dots,n. \]
重心插值公式（对 \(x \neq x_i\)）：
\[p(x) = \frac{\sum_{i=0}^n \frac{w_i}{x - x_i} f_i}{\sum_{i=0}^n \frac{w_i}{x - x_i}}. \]
此形式在数值上更稳定（避免直接计算高次多项式乘积）。为求导数，我们需要对上述有理形式求导。
第四步：推导基于重心权公式的数值微分矩阵
记 \(L(x) = \sum_{j=0}^n \frac{w_j}{x - x_j}\)，则 \(p(x) = \frac{N(x)}{L(x)}\) 其中 \(N(x) = \sum_{i=0}^n \frac{w_i}{x - x_i} f_i\)。
利用商法则求导：
\[p'(x) = \frac{N'(x)L(x) - N(x)L'(x)}{[L(x)]^2}. \]
在节点 \(x = x_k\) 处，由于 \(L(x)\) 和 \(N(x)\) 有奇性，需谨慎处理。可通过极限推导出节点处的导数公式。最终得到重心权微分矩阵元素（对 \(k \neq i\)）：
\[D_{ki} = \ell_i'(x_k) = \frac{w_i / w_k}{x_k - x_i}, \quad k \neq i. \]
对角元素 \(D_{kk}\) 由插值条件 \(\sum_{i=0}^n \ell_i(x) \equiv 1\) 对 \(x\) 求导得：
\[\sum_{i=0}^n \ell_i'(x) = 0 \quad \Rightarrow \quad D_{kk} = -\sum_{\substack{i=0 \\ i \neq k}}^n D_{ki}. \]
因此，只需先计算所有重心权 \(w_i\)（\(O(n^2)\) 次运算，但可预计算），之后 \(D_{ki}\) 的计算仅需 \(O(1)\) 次除法和减法，整个矩阵 \(D\) 可在 \(O(n^2)\) 时间内构造。
第五步：误差分析
数值微分的误差来自插值余项 \(f(x) - p(x)\) 的导数。已知拉格朗日插值余项：
\[f(x) - p(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{j=0}^n (x - x_j). \]
两边求导，代入 \(x = x_k\)，由于右端乘积项在 \(x_k\) 处为零，需用洛必达法则。最终得到在节点 \(x_k\) 处的截断误差：
\[f'(x_k) - p'(x_k) = \frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{\substack{j=0 \\ j \neq k}}^n (x_k - x_j). \]
误差依赖于 \(f^{(n+1)}\) 和节点分布。若节点密集且 \(f\) 光滑，误差随 \(n\) 增大而快速减小。但在高次插值（\(n\) 大）时，可能因龙格现象导致误差振荡，故实际中常用分段低次插值（如分段三次样条）而非全局高次插值。
第六步：在非均匀节点下的高精度实现技巧
- 节点选择优化：对非均匀节点，可选用切比雪夫节点（在区间两端密集）来最小化插值误差，从而提升微分精度。
- 重心权计算的数值稳定性：原始公式 \(w_i = 1 / \prod_{j \neq i}(x_i - x_j)\) 在节点接近时可能上溢/下溢。改用对数求和或缩放技术（如将所有 \(w_i\) 除以同一个常数，不影响 \(D_{ki}\) 比值）提高稳定性。
- 处理端点微分：在区间端点，单侧导数可用非对称的微分矩阵，但精度可能降低。可通过在端点外侧添加虚拟节点（外推）来保持对称性，但需谨慎。
- 高精度计算：使用高精度浮点运算（如四倍精度）可减少舍入误差，尤其当节点非常接近时。
第七步：简单数值示例
设节点 \(x_0=0, x_1=0.5, x_2=1\)，计算 \(f(x)=e^x\) 在 \(x=0.5\) 处的导数近似。
- 计算重心权：
\(w_0 = 1/[(0-0.5)(0-1)] = 1/(0.5) = 2\)
\(w_1 = 1/[(0.5-0)(0.5-1)] = 1/(0.5 \cdot (-0.5)) = -4\)
\(w_2 = 1/[(1-0)(1-0.5)] = 1/(1 \cdot 0.5) = 2\) - 计算 \(D_{1i}\)（在 \(x_1=0.5\) 处）：
\(D_{10} = w_0/(w_1 (x_1-x_0)) = 2/[(-4)(0.5-0)] = 2/(-2) = -1\)
\(D_{12} = w_2/(w_1 (x_1-x_2)) = 2/[(-4)(0.5-1)] = 2/[(-4)(-0.5)] = 2/2 = 1\)
\(D_{11} = - (D_{10}+D_{12}) = -(-1+1) = 0\) - 近似导数：
\(p'(0.5) = (-1) \cdot e^0 + 0 \cdot e^{0.5} + 1 \cdot e^1 = -1 + 0 + 2.71828 = 1.71828\)
精确值 \(f'(0.5)=e^{0.5}\approx 1.64872\)，误差约 0.06956，因节点较少。
总结
本题展示了基于非均匀节点的数值微分的高效构造法：利用重心权公式将微分矩阵表示为简单代数形式，避免了直接对拉格朗日基求导的复杂运算。该方法在节点分布任意时仍保持稳定性，且可通过优化节点分布（如切比雪夫节点）和采用高精度运算进一步提高精度。它适用于实验数据拟合、自适应网格微分等问题，是处理非均匀节点数值微分的有力工具。