C++求一维数组样条差值

一维数组的样条插值可以通过三次样条函数实现。以下是一个简单的 C 代码示例：

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

#define N 10

void spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
{
    double p, qn, sig, un, *u;
    int i, k;

    u = (double *) malloc((n-1) * sizeof(double));

    if (yp1 > 0.99e30) {
        y2[0] = u[0] = 0.0;
    } else {
        y2[0] = -0.5;
        u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0])/(x[1]-x[0]) - yp1);
    }

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

    if (ypn > 0.99e30) {
        qn = un = 0.0;
    } else {
        qn = 0.5;
        un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
    }

    y2[n-1] = (un - qn*u[n-2]) / (qn*y2[n-2] + 1.0);

    for (k = n-2; k >= 0; k--) {
        y2[k] = y2[k]*y2[k+1] + u[k];
    }

    free(u);
}

double splint(double *xa, double *ya, double *y2a, int n, double x)
{
    int klo, khi, k;
    double h, b, a;

    klo = 0;
    khi = n-1;

    while (khi-klo > 1) {
        k = (khi + klo) / 2;

        if (xa[k] > x) {
            khi = k;
        } else {
            klo = k;
        }
    }

    h = xa[khi] - xa[klo];

    if (h == 0.0) {
        fprintf(stderr, "Error: Bad xa input to routine splint\n");
        exit(1);
    }

    a = (xa[khi] - x) / h;
    b = (x - xa[klo]) / h;

    return a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h) / 6.0;
}

int main()
{
    double x[N] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
    double y[N] = {0.1, 0.2, 0.3, 0.5, 0.8, 1.2, 1.7, 2.3, 3.0, 4.0};
    double yp1 = 1e30;
    double ypn = 1e30;
    double y2[N];
    double xi, yi;
    int i;

    spline(x, y, N, yp1, ypn, y2);

    for (i = 0; i < 100; i++) {
        xi = 1.0 + 0.09*i;
        yi = splint(x, y, y2, N, xi);
        printf("%lf %lf\n", xi, yi);
    }

    return 0;
}

这段代码假设您已经有一个包含自变量 x 和因变量 y 的一维数组，以及需要进行插值的点的数量 N。它将在指定的点 xi 处计算样条插值 yi。

`spline()` 函数使用 Thomas 算法计算三次样条函数的二阶导数，而 `splint()` 函数使用这些导数来计算插值点的值。在 `main()` 函数中，我们首先使用 `spline()` 函数计算二阶导数，然后对每个插值点使用 `splint()` 函数计算值，并将结果打印到标准输出。