三次样条插值算法c语言

以下是三次样条插值算法的C语言实现代码：

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

#define N 100

double x[N], y[N], h[N], b[N], u[N], v[N], z[N], c[N], d[N];
int n;

void spline()
{
    int i, k;
    double p, qn, sig, un;

    for (i = 1; i < n; i++)
        h[i] = x[i] - x[i - 1];

    for (i = 2; i < n; i++)
    {
        sig = h[i] / (h[i] + h[i - 1]);
        p = sig * u[i - 1] + 2.0;
        u[i] = (sig - 1.0) / p;
        b[i] = (y[i] - y[i - 1]) / h[i] - (y[i - 1] - y[i - 2]) / h[i - 1];
        b[i] = (6.0 * b[i] / (h[i] + h[i - 1]) - sig * b[i - 1]) / p;
    }

    qn = un = 0.0;
    c[n - 1] = 0.0;

    for (k = n - 2; k >= 1; k--)
    {
        c[k] = u[k] * c[k + 1] + b[k];
        z[k] = (c[k + 1] - c[k]) / h[k];
        qn += (y[k + 1] - y[k]) / (x[k + 1] - x[k]);
        un += z[k] * h[k];
    }

    z[0] = z[1];
    c[0] = 3.0 * qn / h[1] - un * h[1] / 6.0;
    d[0] = -c[0] / h[0];
    d[n - 1] = 0.0;

    for (k = 1; k < n; k++)
    {
        d[k] = (c[k] - c[k - 1]) / h[k];
        b[k] = h[k] * (2.0 * c[k] + c[k - 1]) / 6.0 + (y[k] - y[k - 1]) / h[k];
    }
}

double f(double t)
{
    int i;
    double p;

    i = 0;

    while (x[i + 1] < t)
        i++;

    p = (t - x[i]) / h[i];
    return ((1.0 - p) * y[i] + p * y[i + 1] + p * (1.0 - p) * (c[i] * (1.0 - p) + c[i + 1] * p) * h[i] * h[i] / 6.0);
}

int main()
{
    int i;
    double t;

    printf("请输入数据点个数n：");
    scanf("%d", &n);

    printf("请输入数据点的x坐标和y坐标：\n");

    for (i = 0; i < n; i++)
        scanf("%lf%lf", &x[i], &y[i]);

    spline();

    printf("请输入插值点的x坐标：");
    scanf("%lf", &t);

    printf("在x=%g处的函数值为%g\n", t, f(t));

    return 0;
}