请用C语言写一个光滑不等距插值法

光滑不等距插值法（Smooth non-uniform interpolation）是一种插值方法，用于在非等距离散数据点之间进行插值，以获得连续的函数值。

以下是一个用C语言实现的光滑不等距插值法的代码示例：

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

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

double x[N] = {0.0, 0.5, 1.0, 1.5, 2.0};  // 数据点横坐标
double y[N] = {0.0, 1.2, 1.0, 0.8, 0.0};  // 数据点纵坐标

double phi(double t) {
    return pow(t, 3) * (10.0 + t * (-15.0 + 6.0 * t));
}

double interp(double xi) {
    double h[N], lambda[N], mu[N];
    int i, j, k;
    for (i = 0; i < N; i++) {
        h[i] = (i == 0) ? x[1] - x[0] : x[i] - x[i - 1];
    }
    for (i = 1; i < N - 1; i++) {
        lambda[i] = h[i] / (h[i] + h[i + 1]);
        mu[i] = 1.0 - lambda[i];
    }
    double s[N], t[N], u[N];
    for (i = 1; i < N - 1; i++) {
        s[i] = (y[i] - y[i - 1]) / h[i];
        t[i] = (y[i + 1] - y[i]) / h[i + 1];
        u[i] = 6.0 * (mu[i] * s[i] + lambda[i] * t[i]);
    }
    double A[N], B[N], C[N], D[N], L[N], U[N], Z[N];
    A[0] = 0.0;
    B[0] = 2.0;
    C[0] = 1.0;
    D[0] = 3.0 * (y[1] - y[0]) / h[0] - 3.0 * phi(lambda[1]);
    for (i = 1; i < N - 1; i++) {
        A[i] = mu[i];
        B[i] = 2.0;
        C[i] = lambda[i];
        D[i] = u[i] - mu[i] * u[i - 1];
    }
    A[N - 1] = 1.0;
    B[N - 1] = 2.0;
    C[N - 1] = 0.0;
    D[N - 1] = 3.0 * phi(mu[N - 2]) - 3.0 * (y[N - 1] - y[N - 2]) / h[N - 2];
    for (i = 1; i < N; i++) {
        L[i] = A[i - 1] / B[i - 1];
        U[i] = B[i] - L[i] * C[i - 1];
        Z[i] = (D[i] - L[i] * Z[i - 1]) / U[i];
    }
    double c[N], b[N], d[N];
    c[N - 1] = 0.0;
    for (i = N - 2; i >= 0; i--) {
        c[i] = Z[i + 1] - L[i + 1] * c[i + 1];
        b[i] = (y[i + 1] - y[i]) / h[i] - h[i] * (c[i + 1] + 2.0 * c[i]) / 3.0;
        d[i] = (c[i + 1] - c[i]) / (3.0 * h[i]);
    }
    int l;
    for (l = 0; l < N - 1; l++) {
        if (xi >= x[l] && xi <= x[l + 1]) {
            break;
        }
    }
    double dx = xi - x[l];
    double result = y[l] + b[l] * dx + c[l] * dx * dx + d[l] * dx * dx * dx;
    return result;
}

int main() {
    double xi[M];
    double yi[M];
    double step = (x[N - 1] - x[0]) / (M - 1);
    int i;
    for (i = 0; i < M; i++) {
        xi[i] = x[0] + i * step;
        yi[i] = interp(xi[i]);
    }
    for (i = 0; i < M; i++) {
        printf("%lf %lf\n", xi[i], yi[i]);
    }
    return 0;
}

该代码中的 `interp()` 函数使用了克鲁兹（Krogh）的光滑不等距插值方法，该方法使用三次多项式进行插值，并且要求插值结果在数据点处的导数与原函数相等。为了实现这一要求，我们需要求解一个线性方程组来确定多项式系数。

在 `main()` 函数中，我们生成了 $100$ 个均匀分布在数据点之间的插值点，并使用 `interp()` 函数对这些点进行插值，最终输出插值结果。