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

三次样条插值是一种插值方法，它使用一条光滑的三次函数来逼近一组数据点。以下是一个用C语言编写的三次样条插值程序，其中使用了线性方程组求解器来求解三次样条插值的系数。

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

#define N 10

double x[N], y[N], h[N], b[N], c[N], d[N];

// 解线性方程组 Ax = b，其中 A 是一个三对角矩阵
void solve_tridiagonal(int n, double *a, double *b, double *c, double *d, double *x)
{
    int i;
    double *w = malloc(n * sizeof(double));
    double *v = malloc(n * sizeof(double));
    double beta = b[0];
    x[0] = d[0] / beta;

    for (i = 1; i < n; i++) {
        w[i] = c[i-1] / beta;
        beta = b[i] - a[i] * w[i];
        x[i] = (d[i] - a[i] * x[i-1]) / beta;
    }

    for (i = n-2; i >= 0; i--) {
        v[i] = x[i] - w[i+1] * v[i+1];
    }

    for (i = 0; i < n; i++) {
        x[i] = v[i];
    }

    free(w);
    free(v);
}

void cubic_spline(int n, double *x, double *y)
{
    int i;
    double *a = malloc(n * sizeof(double));
    double *b = malloc(n * sizeof(double));
    double *d = malloc(n * sizeof(double));

    // 计算步长
    for (i = 0; i < n-1; i++) {
        h[i] = x[i+1] - x[i];
    }

    // 计算三对角矩阵的系数
    b[0] = 2 * h[0];
    c[0] = h[0];
    d[0] = 3 * (y[1] - y[0]) / h[0];

    for (i = 1; i < n-1; i++) {
        b[i] = 2 * (h[i-1] + h[i]) - h[i-1] * h[i-1] / b[i-1];
        c[i] = h[i];
        d[i] = 3 * ((y[i+1] - y[i]) / h[i] - (y[i] - y[i-1]) / h[i-1] + h[i-1] * d[i-1] / b[i-1]);
    }

    b[n-1] = 2 * h[n-2] - h[n-2] * h[n-2] / b[n-2];
    d[n-1] = 3 * (y[n-1] - y[n-2]) / h[n-2] - h[n-2] * d[n-2] / b[n-2];

    // 解线性方程组
    solve_tridiagonal(n, a, b, c, d, c);

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

    // 输出插值结果
    printf("x\ty\n");
    for (i = 0; i < n-1; i++) {
        double t;
        for (t = x[i]; t < x[i+1]; t += 0.01) {
            double s = t - x[i];
            double v = y[i] + d[i] * s + c[i] * s * s + a[i] * s * s * s;
            printf("%.2f\t%.2f\n", t, v);
        }
    }

    free(a);
    free(b);
    free(d);
}

int main()
{
    int i, n;
    printf("请输入数据点的个数：");
    scanf("%d", &n);

    printf("请输入数据点的x坐标和y坐标：\n");
    for (i = 0; i < n; i++) {
        scanf("%lf%lf", &x[i], &y[i]);
    }

    cubic_spline(n, x, y);
    return 0;
}

该程序先计算出步长，然后用三对角矩阵的方法求解三次样条插值的系数，最后计算出插值函数的值并输出。