【算法】样条插值算法理论与C++实现

-Aerolite-

已于 2024-11-24 10:39:59 修改

阅读量2.4k

点赞数 22

分类专栏：算法应用文章标签：算法 c++

于 2024-10-26 15:38:25 首次发布

本文链接：https://blog.youkuaiyun.com/weixin_43348948/article/details/143244106

版权

算法应用专栏收录该内容

1 篇文章

订阅专栏

样条插值（Spline Interpolation）是一种在数值分析和计算机图形学中常用的插值方法。通过在已知数据点之间构造一个平滑的曲线（或称为样条），来估计或预测未知数据点的值。样条插值特别适用于需要平滑过渡和连续导数的情况。

一、样条插值的类型

线性样条插值：
- 每个区间内使用线性函数（直线）进行插值。
- 简单，但通常不够平滑。
二次样条插值：
- 每个区间内使用二次多项式（抛物线）进行插值。
- 相比线性样条，可以提供更平滑的曲线，但可能在某些点处导数不连续。
三次样条插值（Cubic Spline Interpolation）：
- 每个区间内使用三次多项式进行插值。
- 提供了平滑的曲线，并且在所有点处导数都是连续的。
- 是最常用的样条插值类型之一。

二、三次样条插值

2.1 关键特性

连续性：插值函数在数据点处是连续的。
平滑性：插值函数的一阶导数（斜率）和二阶导数（曲率）在数据点处也是连续的。
局部性：每个三次多项式只依赖于相邻的三个数据点，因此修改一个数据点只会影响相邻的局部区域。

2.2 实现步骤

三次样条插值（Cubic Spline Interpolation）的核心思想是通过分段三次多项式函数对数据点进行插值。每一段曲线的光滑性和连续性由边界条件和连续性条件决定。

2.2.1 样条函数的定义

对于每个相邻的两点 $x_i, y_i)$ 和 $x_{i+1}, y_{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$ 是需要确定的系数，可以通过已知条件（如函数值、一阶导数和二阶导数在数据点处的连续性）来求解， $S_i(x)$ 在区间 $x_i, x_{i+1}]$ 内有效。

2.2.2 条件与方程

为了确保曲线在每个点连续、光滑，三次样条插值需要满足以下条件：

插值条件：样条曲线必须通过每个已知的数据点：
$S_i(x_i) = y_i \quad \text{和} \quad S_i(x_{i+1}) = y_{i+1}$
这给出了 $a_i = y_i$ 。
连续性条件：曲线在每个内点处的一阶导数和二阶导数必须连续：
- $S'_i(x_{i+1}) = S'_{i+1}(x_{i+1})$ （一阶导数连续）
- $S''_i(x_{i+1}) = S''_{i+1}(x_{i+1})$ （二阶导数连续）
边界条件（自然边界条件）：二阶导数在端点处为零，确保曲线在边界处平滑：
$S''_0(x_0) = 0 \quad \text{和} \quad S''_n(x_n) = 0$

2.2.3 推导系数的过程

（1）计算 $h_i$ 和 $\alpha_i$

定义每个区间的步长 $h_i = x_{i+1} - x_i$ ，并计算中间量 $\alpha_i$ ，它与插值点的值有关：

$\alpha_i = \frac{3}{h_i} (y_{i+1} - y_i) - \frac{3}{h_{i-1}} (y_i - y_{i-1})$

（2）构建三对角矩阵

为了求解未知的 $c_i$ 系数（即每个多项式的二阶导数），我们构建一个线性方程组。这个方程组可以通过求解三对角矩阵来得到。方程组为：
$l_i c_{i-1} + 2 (h_{i-1} + h_i) c_i + h_i c_{i+1} = \alpha_i$
其中：

$l_i$ 和 $h_i$ 是从数据点的步长计算得到的。
$c_0 = c_n = 0$ （自然边界条件）。

（3）解出 $c_i$

通过追赶法等方法，可以解出 $c_i$ 的值。然后使用以下公式求解 $b_i$ 和 $d_i$ ：

$b_i = \frac{y_{i+1} - y_i}{h_i} - \frac{h_i}{3} (2c_i + c_{i+1})$

$d_i = \frac{c_{i+1} - c_i}{3h_i}$

2.2.4 最终插值函数

当所有的系数 $a_i$ 、 $b_i$ 、 $c_i$ 、 $d_i$ 都确定之后，最终的插值函数 $S_i(x)$ 就可以在区间 $x_i, x_{i+1}]$ 内用于插值。

三、三次样条插值的C++实现

3.1 实现步骤

构建三对角矩阵：通过样条方程的连续性和光滑性条件构建三对角矩阵。
解线性方程组：使用追赶法或其他解法求解矩阵。
插值：通过求解出的系数，生成样条插值函数。

3.2 代码

#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

struct Spline {
    vector<double> a, b, c, d, x;
};

// 三次样条插值的构建函数
Spline cubicSpline(const vector<double>& x, const vector<double>& y) {
    int n = x.size() - 1;
    vector<double> a(y), b(n), d(n), h(n), alpha(n), l(n+1), mu(n+1), z(n+1), c(n+1);

    for (int i = 0; i < n; ++i) {
        h[i] = x[i+1] - x[i];
    }

    for (int i = 1; i < n; ++i) {
        alpha[i] = (3.0 / h[i]) * (a[i+1] - a[i]) - (3.0 / h[i-1]) * (a[i] - a[i-1]);
    }

    l[0] = 1.0;
    mu[0] = z[0] = 0.0;

    for (int i = 1; i < n; ++i) {
        l[i] = 2.0 * (x[i+1] - x[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];
    }

    l[n] = 1.0;
    z[n] = c[n] = 0.0;

    for (int j = n-1; j >= 0; --j) {
        c[j] = z[j] - mu[j] * c[j+1];
        b[j] = (a[j+1] - a[j]) / h[j] - h[j] * (c[j+1] + 2.0 * c[j]) / 3.0;
        d[j] = (c[j+1] - c[j]) / (3.0 * h[j]);
    }

    return {a, b, c, d, x};
}

// 使用三次样条插值计算值
double splineEval(const Spline& spline, double x_val) {
    const auto& a = spline.a;
    const auto& b = spline.b;
    const auto& c = spline.c;
    const auto& d = spline.d;
    const auto& x = spline.x;

    int n = x.size() - 1;
    int i = n - 1;

    // 找到插值区间
    for (int j = 0; j < n; ++j) {
        if (x_val >= x[j] && x_val <= x[j+1]) {
            i = j;
            break;
        }
    }

    double dx = x_val - x[i];
    return a[i] + b[i] * dx + c[i] * dx * dx + d[i] * dx * dx * dx;
}

int main() {
    // 样本数据点
    vector<double> x = {0.0, 1.0, 2.0, 3.0, 4.0};
    vector<double> y = {0.0, 1.0, 0.0, 1.0, 0.0};

    // 生成样条插值
    Spline spline = cubicSpline(x, y);

    // 在给定的x值上计算插值
    double x_val = 2.5;
    double y_val = splineEval(spline, x_val);

    cout << "在 x = " << x_val << " 时的插值 y = " << y_val << endl;

    return 0;
}