基于非均匀节点三次样条插值的数值微分公式构造及其截断误差分析
题目描述：
在数值微分中，利用给定的函数在非均匀节点 \(a = x_0 < x_1 < \dots < x_n = b\) 上的值 \(f(x_i)\)，基于三次样条插值构造节点处一阶导数 \(f'(x_i)\) 和二阶导数 \(f''(x_i)\) 的近似计算公式，并分析其截断误差。要求详细阐述三次样条插值在非均匀节点上的构造过程，推导数值微分公式，并进行误差分析。
解题过程循序渐进讲解：
步骤1：问题理解与三次样条插值回顾
数值微分的核心是利用已知离散点近似导数。三次样条插值函数 \(S(x)\) 在每个子区间 \([x_{i-1}, x_i]\) 上是三次多项式，整体具有连续的二阶导数，光滑性好，适合用于导数逼近。在非均匀节点下，节点间距 \(h_i = x_i - x_{i-1}\) 可变，这增加了构造的复杂性。
步骤2：三次样条插值在非均匀节点上的构造
- 定义：在每个子区间 \([x_{i-1}, x_i]\) 上，设样条函数为：
\[S_i(x) = a_i + b_i (x - x_{i-1}) + c_i (x - x_{i-1})^2 + d_i (x - x_{i-1})^3, \quad i=1,\dots,n \]
其中 \(a_i, b_i, c_i, d_i\) 为待定系数。
- 插值条件：\(S_i(x_{i-1}) = f(x_{i-1}), \quad S_i(x_i) = f(x_i)\)。
由此可得：
\[a_i = f(x_{i-1}), \quad a_i + b_i h_i + c_i h_i^2 + d_i h_i^3 = f(x_i) \quad (1) \]
- 连续性条件：在内部节点 \(x_i\) 处，要求 \(S_i'(x_i) = S_{i+1}'(x_i)\) 和 \(S_i''(x_i) = S_{i+1}''(x_i)\)。
记 \(M_i = S''(x_i)\)（称为“矩”），则利用三次多项式的性质，可将系数用 \(M_i\) 表示：
\[c_i = \frac{M_{i-1}}{2}, \quad d_i = \frac{M_i - M_{i-1}}{6h_i} \]
代入(1)解出：
\[b_i = \frac{f(x_i) - f(x_{i-1})}{h_i} - \frac{h_i}{6}(2M_{i-1} + M_i) \]
这样，所有系数由 \(M_i\) 决定。
- 矩 \(M_i\) 的方程组：利用一阶导数连续条件 \(S_i'(x_i) = S_{i+1}'(x_i)\)，整理后得到关于 \(M_i\) 的线性方程组（三对角）：
\[\frac{h_i}{6} M_{i-1} + \frac{h_i + h_{i+1}}{3} M_i + \frac{h_{i+1}}{6} M_{i+1} = \frac{f(x_{i+1}) - f(x_i)}{h_{i+1}} - \frac{f(x_i) - f(x_{i-1})}{h_i} \]
对 \(i=1,\dots,n-1\)。通常补充边界条件（如自然边界 \(M_0 = M_n = 0\) 或固定斜率条件）使方程组封闭。
步骤3：一阶导数近似公式的推导
在节点 \(x_i\) 处，一阶导数 \(f'(x_i)\) 用样条导数 \(S'(x_i)\) 近似。由样条表达式，在区间 \([x_{i-1}, x_i]\) 右端点 \(x_i\) 处有：
\[S_i'(x_i) = b_i + 2c_i h_i + 3d_i h_i^2 \]
代入系数表达式，经代数化简得：
\[f'(x_i) \approx S_i'(x_i) = \frac{f(x_i) - f(x_{i-1})}{h_i} + \frac{h_i}{6}(2M_i + M_{i-1}) \quad (2) \]
类似地，在区间 \([x_i, x_{i+1}]\) 左端点 \(x_i\) 处有对称形式：
\[S_{i+1}'(x_i) = \frac{f(x_{i+1}) - f(x_i)}{h_{i+1}} - \frac{h_{i+1}}{6}(M_{i+1} + 2M_i) \quad (3) \]
在理论上，(2)和(3)应相等（由连续性保证），实际计算时可用其一，或取平均提高精度。
步骤4：二阶导数近似公式的推导
二阶导数直接由矩近似：
\[f''(x_i) \approx M_i \]
因为样条二阶导在节点处连续且正是 \(M_i\)。\(M_i\) 通过解三对角方程组得到，是 \(f''(x_i)\) 的高精度近似。
步骤5：截断误差分析
-
假设：设 \(f \in C^4[a,b]\)，即四阶连续可导。
-
二阶导数的误差：可以证明，对于自然边界条件（\(M_0 = M_n = 0\)），有误差估计：
\[\max_i |f''(x_i) - M_i| = O(h^2) \]
其中 \(h = \max_i h_i\)。推导思路：将差分方程与精确导数 \(f''(x_i)\) 满足的方程相减，利用泰勒展开和矩阵分析得到。
- 一阶导数的误差：以公式(2)为例，将 \(f(x_i), f(x_{i-1})\) 在 \(x_i\) 处泰勒展开，并利用 \(M_{i-1}, M_i\) 与 \(f''(x_{i-1}), f''(x_i)\) 的误差关系，可证得：
\[|f'(x_i) - S_i'(x_i)| = O(h^3) \]
注意，若 \(M_i\) 精确已知（即 \(f''\) 精确），则误差为 \(O(h^3)\)；但因 \(M_i\) 自身有 \(O(h^2)\) 误差，代入后一阶导数整体误差仍为 \(O(h^2)\)。但若节点均匀（\(h_i = h\)），且利用对称性（如(2)和(3)的平均），误差可提升至 \(O(h^4)\)。
- 误差的主要来源：
- 样条插值误差：\(\|f - S\|_\infty = O(h^4)\)
- 导数近似误差：来自样条对导数的逼近阶，通常比函数逼近低一阶。
步骤6：计算流程小结
- 输入非均匀节点 \(x_i\) 及函数值 \(f(x_i)\)。
- 根据边界条件（如自然样条）建立并求解三对角方程组，得到矩 \(M_i\)。
- 用公式(2)或(3)计算一阶导数近似值。
- 直接取 \(M_i\) 作为二阶导数近似值。
- 误差阶：一阶导数 \(O(h^2)\)，二阶导数 \(O(h^2)\)（在非均匀节点下）。
这种方法相比简单的有限差分，在非均匀节点上保持了高阶精度，且数值稳定性更好，尤其适用于导数变化剧烈的函数。实际应用时，可通过节点加密（减小 \(h\) ）进一步提高精度。