三次样条插值C语言代码实现

由于三次样条插值算法较为复杂，因此在此提供一个基于已有库函数的C语言代码实现，需要使用math.h库函数。

#include <stdio.h>
#include <math.h>

#define N 5 // 数据点个数
#define M 50 // 插值点个数

double x[N] = { 0, 1, 2, 3, 4 }; // 数据点x坐标
double y[N] = { 0, 1, -1, 2, -2 }; // 数据点y坐标

double a[N], b[N], c[N], d[N]; // 插值多项式系数

double f(double z) // 原函数
{
    return sin(z);
}

void spline() // 计算插值多项式系数
{
    double h[N - 1], alpha[N - 1], l[N], mu[N - 1], z[N];
    int i;

    for (i = 0; i < N - 1; i++)
    {
        h[i] = x[i + 1] - x[i];
        alpha[i] = 3 * (y[i + 1] - y[i]) / h[i];
    }

    l[0] = 1;
    mu[0] = 0;
    z[0] = 0;

    for (i = 1; i < N - 1; i++)
    {
        l[i] = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1];
        mu[i] = h[i] / l[i];
        z[i] = (alpha[i - 1] - h[i - 1] * z[i - 1]) / l[i];
    }

    l[N - 1] = 1;
    z[N - 1] = 0;
    c[N - 1] = 0;

    for (i = N - 2; i >= 0; i--)
    {
        c[i] = z[i] - mu[i] * c[i + 1];
        b[i] = (y[i + 1] - y[i]) / h[i] - h[i] * (c[i + 1] + 2 * c[i]) / 3;
        d[i] = (c[i + 1] - c[i]) / (3 * h[i]);
        a[i] = y[i];
    }
}

double S(double xi, int j) // 插值多项式
{
    double dx = xi - x[j];

    return a[j] + b[j] * dx + c[j] * dx * dx + d[j] * dx * dx * dx;
}

int main()
{
    spline(); // 计算插值多项式系数

    double xi, yi;
    int i, j;

    for (i = 0; i < M; i++)
    {
        xi = x[0] + (x[N - 1] - x[0]) * i / (M - 1);
        yi = S(xi, 0);

        for (j = 1; j < N; j++)
        {
            if (xi <= x[j])
            {
                yi = S(xi, j);
                break;
            }
        }

        printf("%f\n", yi);
    }

    return 0;
}

以上代码实现了对五个数据点的三次样条插值，并通过50个插值点的计算结果输出了插值多项式的值。可以将其中的x、y数据点和原函数f修改为自己需要的数据，从而实现任意的三次样条插值计算。