基于重心权公式的数值微分及其在非均匀节点下的高精度构造
题目描述
数值微分的目标是近似计算函数 \(f(x)\) 在某点 \(x_0\) 处的导数值 \(f'(x_0)\)。当节点(即已知函数值的点)为非均匀分布时,常见的等距节点差分公式(如三点、五点中心差分)不再直接适用。重心权公式(Barycentric Weights)提供了一种在任意节点分布下高效、稳定地构造插值多项式,并由此推导出高精度数值微分公式的方法。本题要求:
- 理解基于 Lagrange 插值多项式的数值微分公式构造。
- 掌握重心权公式的推导及其优势。
- 针对非均匀节点(如 Chebyshev 节点、Legendre 节点等),利用重心权公式构造数值微分公式,并分析其精度与稳定性。
解题过程
步骤 1:从 Lagrange 插值多项式出发
给定 \(n+1\) 个节点 \(\{x_i\}_{i=0}^n\)(任意分布,不必等距)及函数值 \(f_i = f(x_i)\)。Lagrange 插值多项式为:
\[L_n(x) = \sum_{i=0}^n f_i \, \ell_i(x) \]
其中基函数
\[\ell_i(x) = \prod_{\substack{j=0 \\ j \neq i}}^n \frac{x - x_j}{x_i - x_j} \]
近似导数 \(f'(x) \approx L_n'(x)\)。特别地,在节点 \(x_k\) 处:
\[f'(x_k) \approx \sum_{i=0}^n f_i \, \ell_i'(x_k) \]
定义微分矩阵 \(D_{ki} = \ell_i'(x_k)\),则数值微分转化为矩阵乘法:
\[\mathbf{f}' \approx D \, \mathbf{f} \]
其中 \(\mathbf{f} = [f_0, \dots, f_n]^T\),\(\mathbf{f}'\) 为节点处导数的近似向量。
步骤 2:直接计算微分矩阵的困难
若直接对基函数 \(\ell_i(x)\) 求导,得到公式:
\[\ell_i'(x_k) = \begin{cases} \frac{1}{x_i - x_k} \prod_{\substack{j=0 \\ j \neq i,k}}^n \frac{x_k - x_j}{x_i - x_j}, & i \neq k \\ \sum_{\substack{j=0 \\ j \neq k}}^n \frac{1}{x_k - x_j}, & i = k \end{cases} \]
但此式计算涉及连乘,当 \(n\) 较大时易出现上溢/下溢,且每次计算需 \(O(n^2)\) 操作,效率低。
步骤 3:引入重心权公式简化计算
定义重心权:
\[w_i = \frac{1}{\prod_{\substack{j=0 \\ j \neq i}}^n (x_i - x_j)}, \quad i=0,\dots,n \]
注意 \(w_i\) 与函数值无关,只依赖于节点分布。可预先计算并存储。
利用重心权,Lagrange 基函数可改写为:
\[\ell_i(x) = \frac{ \frac{w_i}{x - x_i} }{ \sum_{j=0}^n \frac{w_j}{x - x_j} } \]
此即重心形式,计算更稳定。
步骤 4:推导基于重心权的微分矩阵公式
对 \(\ell_i(x)\) 求导并整理,可得微分矩阵元素:
\[D_{ki} = \ell_i'(x_k) = \begin{cases} \frac{w_i / w_k}{x_k - x_i}, & i \neq k \\ -\sum_{\substack{j=0 \\ j \neq k}}^n D_{kj}, & i = k \end{cases} \]
说明:
- 当 \(i \neq k\) 时,只需用到预计算的重心权 \(w_i, w_k\) 和节点差 \(x_k - x_i\)。
- 对角元 \(D_{kk}\) 由“零行和”性质(插值常数函数导数为零)确定:\(\sum_{i=0}^n D_{ki} = 0\),故 \(D_{kk} = -\sum_{j \neq k} D_{kj}\)。
优点:
- 避免直接连乘,提高稳定性。
- 计算所有 \(D_{ki}\) 仅需 \(O(n^2)\) 次操作,且常数小。
步骤 5:针对特殊非均匀节点(如 Chebyshev 节点)的重心权
对于 Chebyshev 节点(第一类):
\[x_i = \cos\left( \frac{2i+1}{2n+2} \pi \right), \quad i=0,\dots,n \]
其重心权有显式公式:
\[w_i = (-1)^i \sin\left( \frac{2i+1}{2n+2} \pi \right) \]
(可差一常数因子,因分子分母可约去)。利用此显式可避免连乘计算,进一步保证精度。
步骤 6:精度分析
- 若 \(f(x)\) 是次数 ≤ \(n\) 的多项式,则 Lagrange 插值精确,故数值微分公式在节点处给出精确导数值。
- 对于一般光滑函数,误差来自插值余项。在最大模意义下,Chebyshev 节点使插值误差最小化,因此数值微分精度较高,且通常具有谱精度(即误差随 \(n\) 增加而指数下降,对于解析函数)。
- 误差公式(在节点 \(x_k\) 处):
\[f'(x_k) - \sum_{i=0}^n f_i D_{ki} = \frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{\substack{j=0 \\ j \neq k}}^n (x_k - x_j) \]
步骤 7:稳定性注意事项
- 当 \(n\) 很大时,重心权可能变得极大或极小,但计算 \(D_{ki}\) 时使用的是比值 \(w_i / w_k\),可避免溢出。
- 建议在计算重心权时,先取对数或进行缩放(如令 \(w_0 = 1\) 并相应调整),以控制数值范围。
- 对于等距节点,当 \(n > 20\) 时数值微分往往不稳定(Runge 现象),而非均匀节点(如 Chebyshev)可改善此问题。
步骤 8:算法实现步骤
- 输入节点 \(\{x_i\}_{i=0}^n\) 和函数值 \(\{f_i\}_{i=0}^n\)。
- 计算重心权 \(\{w_i\}\):
- 通用公式:\(w_i = 1 / \prod_{j \neq i} (x_i - x_j)\)。
- 对 Chebyshev 节点,直接用显式 \(w_i = (-1)^i \sin\left( \frac{2i+1}{2n+2} \pi \right)\)。
- 构造微分矩阵 \(D\):
- 对 \(k = 0,\dots,n\):
- 对 \(i \neq k\),计算 \(D_{ki} = (w_i / w_k) / (x_k - x_i)\)。
- 计算 \(D_{kk} = -\sum_{j \neq k} D_{kj}\)。
- 对 \(k = 0,\dots,n\):
- 近似导数向量:\(f'_k \approx \sum_{i=0}^n D_{ki} f_i\)。
步骤 9:示例
用 Chebyshev 节点对 \(f(x) = e^x \sin(5x)\) 在 \([-1,1]\) 上求导:
- 取 \(n=10\) 个节点,计算重心权(显式)。
- 构造微分矩阵 \(D\)。
- 比较近似导数与精确导数 \(f'(x) = e^x (\sin 5x + 5 \cos 5x)\) 在节点处的最大误差。
- 可验证误差随 \(n\) 增加快速衰减(谱精度)。
关键点总结
- 重心权公式将 Lagrange 基的导数计算转化为简单的权重比值,避免了直接连乘的不稳定性。
- 对非均匀节点(尤其是 Chebyshev 节点),重心权有显式表达式,可高效、高精度地构造微分矩阵。
- 该方法尤其适用于高精度要求且节点可自由选择的场景,如谱方法中的配点求导。
- 注意数值稳定性:预先缩放重心权,并避免大 \(n\) 时矩阵 \(D\) 的条件数恶化(可通过选择节点如 Chebyshev 或 Legendre 来缓解)。
通过以上步骤,你应能理解如何利用重心权公式在非均匀节点下构造高精度数值微分公式,并应用于实际计算。