### 手动实现三次样条插值 为了生成三次样条插值曲线并绘制其图形，可以按照以下方法手动计算和实现。以下是详细的理论基础以及代码实现。 #### 节点定义与基函数构造 对于给定的节点向量 \( \mathbf{t} = \{0, 0, 0, 1, 2, 3, 4, 4, 4\} \)[^1]，可以通过递归方式构建 B 样条基函数。具体来说，\( k=4 \) 表示我们需要构造四次（即三次多项式）B 样条基函数。这些基函数通过分段线性和递推关系来定义： \[ N_{i,k}(u) = \frac{u - t_i}{t_{i+k-1} - t_i}N_{i,k-1}(u) + \frac{t_{i+k}-u}{t_{i+k}-t_{i+1}}N_{i+1,k-1}(u), \] 其中边界条件为当 \( u \in [t_i, t_{i+1}) \)，则 \( N_{i,1}(u) = 1 \); 否则为零[^1]。 --- #### 数据准备 已知目标函数为 \( f(x) = \frac{1}{1+x^2} \)，区间为 [-5, 5] 的均匀分布节点数为 11，则对应的节点坐标为： \[ x = \{-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5\}, \] 对应的目标函数值为： \[ y = \left\{\frac{1}{1+(-5)^2}, \dots, \frac{1}{1+(5)^2}\right\}. \] --- #### 构造三对角矩阵求解二阶导数值 三次样条插值的核心在于解决自然样条方程组。设 \( S''(x_i) = M_i \)，那么满足如下差分方程： \[ h_{i-1}M_{i-1} + 2(h_{i-1} + h_i)M_i + h_iM_{i+1} = 6(f[x_{i-1},x_i,x_{i+1}] - f[x_{i-1},x_i]), \] 其中 \( h_i = x_{i+1} - x_i \)，且 \( f[x_i, x_j] \) 是两点之间的差商。 最终得到一个关于 \( M_1, M_2, ..., M_n \) 的三对角线性方程组，可利用 Thomas 算法高效求解。 --- #### Python 实现代码 下面是完整的 Python 实现代码，用于生成三次样条插值曲线及其绘图。 ```python import numpy as np import matplotlib.pyplot as plt def cubic_spline_interpolation(x_data, y_data): n = len(x_data) h = np.diff(x_data) # 差商计算 alpha = np.zeros(n - 1) for i in range(1, n - 1): alpha[i] = (3 / h[i]) * (y_data[i + 1] - y_data[i]) - \ (3 / h[i - 1]) * (y_data[i] - y_data[i - 1]) # 初始化三对角矩阵系数 l = np.ones(n) mu = np.zeros(n) z = np.zeros(n) for i in range(1, n - 1): l[i] = 2 * (x_data[i + 1] - x_data[i - 1]) - h[i - 1] * mu[i - 1] mu[i] = h[i] / l[i] z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i] b = np.zeros(n) c = np.zeros(n) d = np.zeros(n) for j in range(n - 2, -1, -1): # Back substitution c[j] = z[j] - mu[j] * c[j + 1] b[j] = (y_data[j + 1] - y_data[j]) / h[j] - h[j] * (c[j + 1] + 2 * c[j]) / 3 d[j] = (c[j + 1] - c[j]) / (3 * h[j]) return b[:-1], c[:], d[:-1] # 定义输入数据点 x_points = np.linspace(-5, 5, 11) y_points = 1 / (1 + x_points ** 2) b_coeffs, c_coeffs, d_coeffs = cubic_spline_interpolation(x_points, y_points) # 插值函数评估 def spline_eval(x_new, x_data, y_data, b, c, d): result = [] for xi in x_new: idx = np.searchsorted(x_data, xi) - 1 if idx >= len(x_data) - 1 or idx < 0: continue hi = x_data[idx + 1] - x_data[idx] yi = y_data[idx] si = yi + b[idx]*(xi-x_data[idx]) + c[idx]*((xi-x_data[idx])**2)/2 + d[idx]*((xi-x_data[idx])**3)/(6*hi) result.append(si) return np.array(result) # 绘制图像 x_fine = np.linspace(min(x_points), max(x_points), 500) y_fine = spline_eval(x_fine, x_points, y_points, b_coeffs, c_coeffs, d_coeffs) plt.figure(figsize=(8, 6)) plt.plot(x_fine, y_fine, label="Cubic Spline Interpolation", color='blue') plt.scatter(x_points, y_points, color='red', label="Data Points") plt.title("Cubic Spline Interpolation of $f(x) = \\frac{1}{1+x^2}$") plt.xlabel("$x$") plt.ylabel("$f(x)$") plt.legend() plt.grid(True) plt.show() ``` --- #### 结果分析 上述程序实现了基于手动生成的三次样条插值算法，并将其应用于指定的数据集上。可以看到，在较平滑区域（如靠近原点附近），插值效果较好；而在端点处由于高曲率变化较大，可能会存在一定的误差。这种现象可通过调整边界条件进一步优化。 ---