基于python的三次样条插值原理及代码

news2024/12/28 17:50:15

1 三次样条插值

1.1 三次样条插值的基本概念

三次样条插值是通过求解三弯矩方程组（即三次样条方程组的特殊形式）来得出曲线函数组的过程。在实际计算中，还需要引入边界条件来完成计算。样条插值的名称来源于早期工程师制图时使用的细长木条（样条），这些木条被固定在样点上，然后自由弯曲以绘制出曲线。

1.2 三次样条插值的条件

假设在区间[a, b]上有n+1个样本点，即[a, b]区间被划分成n个小区间 $[(x_{0},x_{1}),(x_{1}, x_{2}), ...,(x_{n-1}, x_{n})]$ ，其中 $x_{0}=a$ , $x_{n}=b$ 。三次样条插值需要满足以下条件：

在每个分段区间 $[x_{i},x_{i+1}]$ 上，插值函数 $S(x) = S_{i}(x)$ 都是一个三次多项式。
满足插值条件，即 $S(x_i) = y_i \quad (i = 0, 1, ..., n)$ 。
曲线光滑，即S(x)、S'(x)、S''(x)连续。

1.3 三次样条插值的公式推导

1. 3.1 三次样条函数的形式

在每个小区间 $[x_{i},x_{i+1}]$ 上，三次样条函数 $S_i(x)$ 可以表示为：

$S_i(x) = a_i + b_i(x - x_i) + c_i(x - x_i)2 + d_i(x - x_i)3$

其中， $a_i ,b_i,c_i,d_i$ 是待求的系数。

1.3.2 插值条件

由于 $S(x_{i}) = y_{i}$ ，我们可以得到n+1个方程：

$S_i(x_i) = a_i = y_i \quad (i = 0, 1, ..., n)$

1.3.3 连续性条件

除了插值条件外，还需要保证曲线在样本点处的一阶导数和二阶导数连续。即：

一阶导数连续：对于任意i（0 < i < n），有

$S_i'(x_{i+1}) = S_{i+1}'(x_{i+1})$
二阶导数连续：对于任意i（0 < i < n），有

$S_i''(x_{i+1}) = S_{i+1}''(x_{i+1})$

1.3.4 求解未知数

每个小区间内有4个未知数（ $a_i ,b_i,c_i,d_i$ ），共有n个小区间，所以一共有4n个未知数。根据插值条件、一阶导数连续和二阶导数连续条件，我们可以列出以下方程组：

插值条件：提供n+1个方程（其中ai=yi已知）。
一阶导数连续：提供n-1个方程（两个端点不参与此条件）。
二阶导数连续：提供n-1个方程（同样，两个端点不参与此条件）。

这样，我们总共有4n-2个方程，但还需要两个额外的方程来求解所有未知数。这两个额外的方程由边界条件给出。

1.3.5 边界条件

常用的边界条件有三种：

自然边界（Natural Spline）：指定端点二阶导数为0，即

$S''(x_0) = 0 = S''(x_n)$
固定边界（Clamped Spline）：指定端点一阶导数，即

$S'(x_0) = A, \quad S'(x_n) = B$
非扭结边界（Not-a-Knot Spline）：使两端点的三阶导与这两端点的邻近点的三阶导相等，即

$S'''(x_0) = S'''(x_1), \quad S'''(x_{n-1}) = S'''(x_n)$

1.3.6 求解过程

以自然边界条件为例，即 $S''(x_0) = 0, \quad S''(x_n) = 0$ ，我们需要通过这些条件来求解每个小区间上的二阶导数 $m_{i}=S''(x_i)$ （通常称为样条曲线的“弯矩”或“斜率”）。

步骤一：设置方程

在每个小区间 $[x_{i},x_{i+1}]$ 上，三次样条 $S_i(x)$ 的二阶导数 $S''(x_n)$ 是一个常数 $m_{i}$ 。因此，三次样条 $S_i(x)$ 可以表示为：

$S_i(x) = \frac{m_i}{6}(x_{i+1} - x)3 + \frac{m_{i+1}}{6}(x - x_i)3 + a_i(x_{i+1} - x) + b_i(x - x_i) + c_i$

其中 $a_i ,b_i,c_i$ 是需要通过插值条件和其他连续性条件求解的系数。

步骤二：利用插值条件

由于 $S_i(x_{i}) = y_{i},S_i(x_{i+1}) = y_{i+1}$ ，我们可以列出 n+1 个插值方程。但由于 $a_{i}=b_{i}$ 由 $S_i(x_{i}) = y_{i}$ 直接得出），我们实际上只需要解出 $a_{i},b_{i}$ （或者更直接地，解出 $m_{i}$ ）。

步骤三：利用一阶导数连续性

在 $x=x_{i}+1$ 处，有 $S_i{}'(x_{i}) = S_{i+1}{}'(x_{i+1})$ 。对 $S_i{}(x_{i}) = S_{i+1}{}(x)$ 求导并设置相等，我们可以得到一个关于 $m_{i}$ 和 $m_{i+1}$ 的方程。由于有 n−1 个这样的内部节点，因此我们可以得到 n−1 个这样的方程。

步骤四：利用二阶导数连续性（自然边界条件）

在 $x=x_{0}$ 和 $x=x_{n}$ 处，有 $S''(x_0) = 0, \quad S''(x_n) = 0$ ，即 $m_{0}=0$ 和 $m_{n}=0$ 这两个条件提供了额外的两个方程。

步骤五：解线性方程组

现在，我们有了 n−1 个一阶导数连续性方程，加上两个自然边界条件方程，共 n+1 个方程，用于求解 n 个未知数 $m_{1},m_{2},...m_{n-1}$ （注意 $m_{0}$ 和 $m_{n}$ 已知为0）。这是一个线性方程组，可以通过高斯消元法、LU分解等方法求解。

步骤六：计算系数 $a_i ,b_i,c_i$

一旦我们得到了所有的 mi，就可以通过插值条件和其他连续性条件回代求解出每个小区间上的 $a_i ,b_i,c_i$ 。具体来说，可以使用以下公式：

$a_{i}=b_{i}$ （已知）
$b_{i},c_{i}$ 可以通过 $S_i(x_{i+1}) = y_{i+1}$ 和 $S_i{}'(x_{i}) = S_{i+1}{}'(x_{i+1})$ 的联立方程求解。

步骤七：构建完整的三次样条函数

最后，我们得到了每个小区间上的三次多项式 $S_i(x)$ ，它们共同构成了完整的三次样条插值函数 S(x)。在任意点 x∈ $[x_{i},x_{i+1}]$ 上， $S(x) = S_{i}(x)$ 。

1.4 总结

三次样条插值通过求解一系列线性方程组来确定每个小区间上的三次多项式，这些多项式在样本点处满足插值条件，并且具有连续的一阶导数和二阶导数。自然边界条件、固定边界条件或非扭结边界条件等不同的边界条件选择会影响求解过程和最终的插值曲线形状。三次样条插值因其光滑性和易于实现的特性，在数值分析、数据可视化、计算机辅助设计等领域有着广泛的应用。

2 代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline

# 定义数据点
x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15])
y1 = np.array([0, 4, 6, 8, 9, 11, 13, 16, 19, 21, 23, 25, 28, 30, 33, 34])
y2 = np.array([0, 1, 2, 4, 5, 6, 7, 8, 10, 12, 14, 16, 20, 25, 29, 34])
y3 = np.array([0, 8, 12, 18, 20, 26, 30, 32, 33, 34, 34, 36, 36, 35, 35, 34])

# 创建三次样条插值函数
cs1 = CubicSpline(x, y1)
cs2 = CubicSpline(x, y2)
cs3 = CubicSpline(x, y3)

# 生成插值后的x坐标
x_new = np.linspace(0, 15, 100)

# 计算插值后的y坐标
y_new1 = cs1(x_new)
y_new2 = cs2(x_new)
y_new3 = cs3(x_new)

# 设置颜色选择器和颜色映射
select2 = (0.1, 0.4, 0.5, 0.3, 0.8, 1.0) # 非连续型色组在0-(色组长度-1)之间选择
colors = plt.get_cmap('rainbow')(select2) # 从色组里选择颜色

# 绘制原始数据点和插值曲线
# 调整图表大小和位置
plt.figure(figsize=(10, 5))
plt.plot(x, y1, 'o', color=colors[0], label='ori_y1')
plt.plot(x, y2, 'o', color=colors[2], label='ori_y2')
plt.plot(x, y3, 'o', color=colors[4], label='ori_y3')
plt.plot(x_new, y_new1, '-', color=colors[3], label='y1_interp_curve')
plt.plot(x_new, y_new2, '-', color=colors[1], label='y2_interp_curve')
plt.plot(x_new, y_new3, '-', color=colors[5], label='y3_interp_curve')

# 设置图例和标题
plt.legend()
plt.xticks(np.arange(0, 16, 1))
plt.title('Cubic Spline Interpolation Example')
plt.xlabel('x')
plt.ylabel('y')


# 显示图表
plt.show()