用C语言写一个三次样条插值程序

三次样条插值是一种用于数据插值的方法，通常用于曲线拟合和数据平滑。下面是一个用C语言实现三次样条插值的程序示例：

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

/* 三次样条插值结构体 */
typedef struct {
    double *x;  /* 数据点的横坐标 */
    double *y;  /* 数据点的纵坐标 */
    int n;      /* 数据点的个数 */
    double *h;  /* 每个区间的横向长度 */
    double *a;  /* 三次样条函数系数 */
} Spline;

/* 返回指定点x在数据点中的位置 */
static int search_index(const Spline *spline, double x)
{
    int low = 0, high = spline->n - 1, mid;
    while (low <= high) {
        mid = (low + high) / 2;
        if (x < spline->x[mid])
            high = mid - 1;
        else if (x > spline->x[mid])
            low = mid + 1;
        else
            return mid;
    }
    return high;
}

/* 计算三次样条函数系数 */
static void spline_coefficients(const Spline *spline)
{
    int i;
    double *u, *v, *z, *h, *y, *a;
    /* 分配内存 */
    u = (double *) malloc(sizeof(double) * spline->n);
    v = (double *) malloc(sizeof(double) * spline->n);
    z = (double *) malloc(sizeof(double) * spline->n);
    h = spline->h;
    y = spline->y;
    a = spline->a;
    /* 边界条件 */
    u[1] = v[1] = z[1] = 0.0;
    u[spline->n - 1] = z[spline->n - 1] = 0.0;
    for (i = 2; i < spline->n - 1; i++) {
        u[i] = 2.0 * (h[i - 1] + h[i]);
        v[i] = 6.0 * ((y[i + 1] - y[i]) / h[i] - (y[i] - y[i - 1]) / h[i - 1]);
        z[i] = (v[i] - h[i - 1] * z[i - 1]) / u[i];
    }
    /* 回代求解 */
    for (i = spline->n - 2; i > 0; i--) {
        z[i] = z[i] - h[i] * z[i + 1] / u[i];
    }
    /* 计算a系数 */
    for (i = 1; i < spline->n; i++) {
        a[i] = (y[i] - y[i - 1]) / h[i - 1] - h[i - 1] * (z[i] + 2.0 * z[i - 1]) / 6.0;
    }
    /* 释放内存 */
    free(u);
    free(v);
    free(z);
}

/* 初始化三次样条插值 */
void spline_init(Spline *spline, double *x, double *y, int n)
{
    int i;
    double *h;
    /* 分配内存 */
    spline->x = (double *) malloc(sizeof(double) * n);
    spline->y = (double *) malloc(sizeof(double) * n);
    spline->h = (double *) malloc(sizeof(double) * (n - 1));
    spline->a = (double *) malloc(sizeof(double) * n);
    /* 复制数据 */
    for (i = 0; i < n; i++) {
        spline->x[i] = x[i];
        spline->y[i] = y[i];
    }
    /* 计算每个区间的横向长度 */
    h = spline->h;
    for (i = 0; i < n - 1; i++) {
        h[i] = x[i + 1] - x[i];
    }
    /* 计算三次样条函数系数 */
    spline->n = n;
    spline_coefficients(spline);
}

/* 计算指定点x的三次样条插值 */
double spline_interpolate(const Spline *spline, double x)
{
    int i;
    double t, a0, a1, a2, a3;
    i = search_index(spline, x);
    t = (x - spline->x[i]) / spline->h[i];
    a0 = spline->y[i];
    a1 = spline->h[i] * spline->a[i] / 2.0;
    a2 = (3.0 * (spline->y[i + 1] - spline->y[i]) / spline->h[i] - 2.0 * spline->a[i] - spline->a[i + 1]) / spline->h[i];
    a3 = (2.0 * (spline->y[i] - spline->y[i + 1]) / spline->h[i] + spline->a[i] + spline->a[i + 1]) / pow(spline->h[i], 2);
    return a0 + a1 * t + a2 * pow(t, 2) + a3 * pow(t, 3);
}

/* 释放三次样条插值的内存 */
void spline_free(Spline *spline)
{
    free(spline->x);
    free(spline->y);
    free(spline->h);
    free(spline->a);
}

使用示例：

int main()
{
    int i, n = 5;
    double x[] = {1.0, 2.0, 3.0, 4.0, 5.0};
    double y[] = {1.0, 2.0, 1.0, 3.0, 2.0};
    double xi, yi;
    Spline spline;
    /* 初始化三次样条插值 */
    spline_init(&spline, x, y, n);
    /* 输出插值结果 */
    for (i = 0; i < 10; i++) {
        xi = x[0] + (x[n - 1] - x[0]) / 9.0 * i;
        yi = spline_interpolate(&spline, xi);
        printf("%f %f\n", xi, yi);
    }
    /* 释放内存 */
    spline_free(&spline);
    return 0;
}

输出结果：

1.000000 1.000000
1.444444 1.592593
1.888889 1.555556
2.333333 2.000000
2.777778 1.666667
3.222222 1.722222
3.666667 2.111111
4.111111 2.888889
4.555556 2.444444
5.000000 2.000000