编辑代码

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

#define N 100  // 定义点的数量

// 定义一个结构体来存储每个点的坐标
typedef struct {
    double x;  // x坐标
    double y;  // y坐标
} Point;

// 三次样条插值函数
void cubicSpline(Point points[], int n, double a[], double b[], double c[], double d[]) {
    double h[N - 1], alpha[N - 1], l[N], mu[N], z[N];

    // 计算每个区间的宽度h和斜率alpha
    for (int i = 0; i < n - 1; i++) {
        h[i] = points[i + 1].x - points[i].x;  // 计算区间宽度
        alpha[i] = (points[i + 1].y - points[i].y) / h[i];  // 计算斜率
    }

    l[0] = 1;  // 设置边界条件
    mu[0] = 0; // 设置边界条件
    z[0] = 0;  // 设置边界条件

    // 前向消元过程
    for (int i = 1; i < n - 1; i++) {
        l[i] = 2 * (points[i + 1].x - points[i - 1].x) - h[i - 1] * mu[i - 1];  // 计算l[i]
        mu[i] = h[i] / l[i];  // 计算mu[i]
        z[i] = (alpha[i] - alpha[i - 1] * h[i - 1]) / l[i];  // 计算z[i]
    }

    l[n - 1] = 1;  // 设置边界条件
    z[n - 1] = 0;  // 设置边界条件
    c[n - 1] = 0;  // 设置边界条件

    // 后向替换过程
    for (int j = n - 2; j >= 0; j--) {
        c[j] = z[j] - mu[j] * c[j + 1];  // 计算c[j]
        b[j] = alpha[j] - h[j] * (c[j + 1] + 2 * c[j]) / 3;  // 计算b[j]
        d[j] = (c[j + 1] - c[j]) / h[j];  // 计算d[j]
    }
}

// 打印样条系数的函数
void printSplineCoefficients(int n, double a[], double b[], double c[], double d[]) {
    for (int i = 0; i < n - 1; i++) {
        printf("Segment %d: a = %lf, b = %lf, c = %lf, d = %lf\n", i, a[i], b[i], c[i], d[i]);
    }
}

double cubicSplineEval(Point points[], double b[], double c[], double d[], int n, double xe) {
    int i;
    // 找到xe所在的区间[x[i], x[i+1]]  
    for (i = 0; i < n - 1; i++) {
        if (xe >= points[i].x && xe < points[i+1].x) {
            break;
        }
    }
    // 计算插值点的函数值  
    double h = points[i + 1].x - points[i].x;
    double a = points[i].y;
    double xe_minus_xi = xe - points[i].x;
    return a + b[i] * xe_minus_xi + c[i] * xe_minus_xi * xe_minus_xi + d[i] * xe_minus_xi * xe_minus_xi * xe_minus_xi;
}



int main() {
    Point points[N];  // 定义一个点数组来存储100个点
    double a[N], b[N], c[N], d[N];  // 定义数组来存储样条的系数

    // 生成100个点（这里使用y = x^2作为示例）
    for (int i = 0; i < N; i++) {
        points[i].x = (double)i;  // x从0到10
        points[i].y = (points[i].x * points[i].x);  // y = x^2
    }

    // 计算三次样条的系数
    cubicSpline(points, N, a, b, c, d);

    // 打印系数
    printSplineCoefficients(N, a, b, c, d);

    printf("Cubic Spline Interpolation results:\n");  
    for (int i = 0; i < N; i += 1) {  
       printf("f(%5.2f) = %8.4f %8.4f %8.4f\n", points[i].x+0.5, cubicSplineEval(points, b, c, d, N, points[i].x+0.5), ((points[i].x + 0.5)* (points[i].x + 0.5)), (cubicSplineEval(points, b, c, d, N, points[i].x + 0.5)- ((points[i].x + 0.5) * (points[i].x + 0.5))));
     } 

    return 0;  // 程序结束
}