用C++实现一份三次样条插值

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

#define N 10

double x[N + 1] = {0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0};
double y[N + 1] = {1.0, 1.6, 2.0, 2.2, 2.2, 2.0, 1.6, 1.0, 0.2, -0.6, -1.0};
double m[N + 1];

double h(int i) {
return x[i + 1] - x[i];
}

double A(int i) {
return h(i) / 6.0;
}

double B(int i) {
return (x[i + 1] + x[i]) / 2.0;
}

double C(int i) {
return h(i) / 3.0;
}

double D(int i) {
return (y[i + 1] - y[i]) / h(i) - m[i + 1] * h(i) / 6.0 - m[i] * h(i) / 3.0;
}

void init_m() {
double* a = (double*)malloc(sizeof(double) * (N - 1));
double* b = (double*)malloc(sizeof(double) * (N - 1));
double* c = (double*)malloc(sizeof(double) * (N - 1));
double* d = (double*)malloc(sizeof(double) * (N - 1));
double* l = (double*)malloc(sizeof(double) * (N + 1));
double* u = (double*)malloc(sizeof(double) * (N + 1));
double* z = (double*)malloc(sizeof(double) * (N + 1));

for (int i = 1; i < N; i++) {
a[i - 1] = A(i - 1);
b[i - 1] = -A(i - 1) - A(i);
c[i - 1] = A(i);
d[i - 1] = D(i - 1);
}

l[0] = 1.0;
u[0] = 0.0;
z[0] = 0.0;

for (int i = 1; i < N; i++) {
l[i] = 2.0 * (x[i + 1] - x[i - 1]) - h(i - 1) * u[i - 1];
u[i] = h(i) / l[i];
z[i] = (d[i - 1] - h(i - 1) * z[i - 1]) / l[i];
}

l[N] = 1.0;
z[N] = 0.0;
m[N] = 0.0;

for (int i = N - 1; i >= 0; i--) {
m[i] = z[i] - u[i] * m[i + 1];
}

free(a);
free(b);
free(c);
free(d);
free(l);
free(u);
free(z);
}

double S(int i, double t) {
double A = (m[i + 1] - m[i]) / (6.0 * h(i));
double B = m[i] / 2.0;
double C = (y[i + 1] - y[i]) / h(i) - h(i) * (2.0 * m[i] + m[i + 1]) / 6.0;
double D = y[i] - m[i] * h(i) * h(i) / 6.0;

return A * pow(t - x[i], 3) + B * pow(t - x[i], 2) + C * (t - x[i]) + D;
}

double f(double t) {
int i = 0;
while (i < N && x[i + 1] < t) {
i++;
}
if (i == N) {
i--;
}
return S(i, t);
}

int main() {
init_m();

printf("x\t\ty\t\tf(x)\n");
for (double t = 0.0; t <= 2.0; t += 0.1) {
printf("%lf\t%lf\t%lf\n", t, f(t), y[(int)(t * 5)]);
}

return 0;
}