三次样条插值算法曲线拟合c代码

三次样条插值算法是一种常用的曲线拟合算法。它通过将给定数据点之间的曲线分段拟合为三次函数，从而得到一个平滑的曲线拟合结果。下面给出一个三次样条插值算法的c代码实现。

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

#define MAX_POINTS 100 // 最大数据点数目

// 数据点结构体
typedef struct _point {
double x;
double y;
} point;

// 三次函数结构体
typedef struct _cubic {
double a, b, c, d;
} cubic;

static point data[MAX_POINTS]; // 数据点数组
static int n = 0; // 数据点数目

// 初始化数据点
void init_data() {
data[n++] = (point) { 0.0, 0.0 };
data[n++] = (point) { 1.0, 1.0 };
data[n++] = (point) { 2.0, 0.0 };
}

// 求解三次函数系数
void solve_cubic(cubic *c, double x0, double x1, double y0, double y1, double dy0, double dy1) {
double dx = x1 - x0;
c->a = (dy0 + dy1 - 2*(y0 - y1)/dx) / (dx*dx);
c->b = (3*(y0 - y1)/dx - 2*dy0 - dy1) / dx;
c->c = dy0;
c->d = y0;
}

// 计算三次函数的值
double eval_cubic(cubic *c, double x) {
double dx = x - c->d;
return c->a*dx*dx*dx + c->b*dx*dx + c->c*dx + c->d;
}

// 三次样条插值算法曲线拟合
void cubic_spline() {
int i;
cubic *c = (cubic*)malloc((n-1) * sizeof(cubic));
double *h = (double*)malloc(n * sizeof(double));
double *alpha = (double*)malloc(n * sizeof(double));
double *l = (double*)malloc((n-1) * sizeof(double));
double *u = (double*)malloc((n-1) * sizeof(double));
double *z = (double*)malloc(n * sizeof(double));
double *b = (double*)malloc((n-1) * sizeof(double));
double *c_ = (double*)malloc(n * sizeof(double));
double *d = (double*)malloc((n-1) * sizeof(double));

// 计算 h_i
for (i = 0; i < n-1; i++)
h[i] = data[i+1].x - data[i].x;

// 计算 alpha_i
for (i = 1; i < n-1; i++)
alpha[i] = 3/h[i]*(data[i+1].y - data[i].y) - 3/h[i-1]*(data[i].y - data[i-1].y);

// 计算 l,u,z
l[0] = 1;
u[0] = z[0] = 0;
for (i = 1; i < n-1; i++) {
l[i] = 2*(data[i+1].x - data[i-1].x) - h[i-1]*u[i-1];
u[i] = h[i]/l[i];
z[i] = (alpha[i] - h[i-1]*z[i-1])/l[i];
}
l[n-1] = 1;
z[n-1] = c_[n-1] = 0;
for (i = n-2; i >= 0; i--) {
c_[i] = z[i] - u[i]*c_[i+1];
b[i] = (data[i+1].y - data[i].y)/h[i] - h[i]*(c_[i+1] + 2*c_[i])/3;
d[i] = (c_[i+1] - c_[i])/(3*h[i]);
solve_cubic(&c[i], data[i].x, data[i+1].x, data[i].y, data[i+1].y, b[i], b[i+1]);
}

// 输出拟合结果
printf("拟合结果：\n");
for (i = 0; i < n-1; i++) {
printf("[%.2f,%.2f], (%.2f*x^3 + %.2f*x^2 + %.2f*x + %.2f)\n",
data[i].x, data[i+1].x, c[i].a, c[i].b, c[i].c, c[i].d);
}
printf("\n");

free(c);
free(h);
free(alpha);
free(l);
free(u);
free(z);
free(b);
free(c_);
free(d);
}

int main() {
init_data();
cubic_spline();
return 0;
}