编辑代码

#include <stdio.h>  
#include <stdlib.h>  
  
#define N 100  
  
void spline(double x[], double y[], double b[], double c[], double d[], int n) {  
    double h;  
    int i;  
    double *alpha = malloc((n - 1) * sizeof(double));  
    double *l = malloc(n * sizeof(double));  
    double *mu = malloc(n * sizeof(double));  
    double *z = malloc(n * sizeof(double));  
  
    // Step 1: Compute the divided differences (slopes) and set up the tridiagonal system  
    for (i = 0; i < n - 1; i++) {  
        h = x[i + 1] - x[i];  
        alpha[i] = (3 / h) * (y[i + 1] - y[i]) - (3 / h) * (y[i] - y[i - 1]);  
    }  
  
    // Set up the boundary conditions (natural spline)  
    l[0] = 1;  
    mu[0] = 0;  
    z[0] = 0;  
  
    // Step 2: Solve the tridiagonal system using the LU decomposition method  
    for (i = 1; i < n - 1; i++) {  
        h = x[i + 1] - x[i];  
        l[i] = 2 * (x[i + 1] - x[i - 1]) - h * mu[i - 1];  
        mu[i] = h / l[i];  
        z[i] = (alpha[i] - h * z[i - 1]) / l[i];  
    }  
  
    l[n - 1] = 1; // Boundary condition for the last point  
  
    // Step 3: Compute the coefficients b[i], c[i], d[i] from the LU solution  
    for (i = n - 2; i >= 0; i--) {  
        c[i] = z[i] - mu[i] * c[i + 1];  
        b[i] = alpha[i] - (x[i + 1] - x[i]) * (c[i + 1] + 2 * c[i]) / 3;  
        d[i] = (c[i + 1] - c[i]) / (3 * (x[i + 1] - x[i]));  
    }  
  
    // Set the last coefficients to zero (natural spline)  
    b[n - 1] = 0;  
    c[n - 1] = 0;  
    d[n - 1] = 0;  
  
    free(alpha);     // 释放动态分配的内存  
    free(l);  
    free(mu);  
    free(z);  
}  
  
// 用于计算三次样条插值在给定点的值的函数  
double cubicSplineEval(double x[], double y[], 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 >= x[i] && xe < x[i + 1]) {  
            break;  
        }  
    }  
    // 计算插值点的函数值  
    double h = x[i + 1] - x[i];  
    double a = y[i];  
    double xe_minus_xi = xe - x[i];  
    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() {  
    // 示例数据  
    double x[N] ;  
    double y[N] ;  
    double b[N], c[N], d[N]; // 存储三次样条插值的系数  
  
    // 初始化数据点  
    for (int i = 0; i < N; i++) {  
        x[i] = i;  
        y[i] = x[i] * x[i] - x[i] + 3;  // 示例函数：f(x) = x^2 - x + 3  
    }  
    // 计算三次样条插值的系数  
    spline(x, y, b, c, d, N);  
    
    // 输出插值结果  
//    printf("Cubic Spline Interpolation results:\n");  
//    for (int i = 0; i < N; i += 10) {  
//       printf("f(%5.2f) = %8.4f\n", x[i], cubicSplineEval(x, y, b, c, d, N, x[i]));  
//   } 
    double xe; // 要计算插值的x值  
    // 假设我们要计算xe=某个值时的插值结果  
    xe = 81.0;/* 填入要计算的x值 */;  
    double ye = cubicSplineEval(x, y, b, c, d, N, xe);  
    printf("在x=%f处的插值结果是: %f\n", xe, ye);  

    xe = 81.6;/* 填入要计算的x值 */;  
    ye = cubicSplineEval(x, y, b, c, d, N, xe);  
    printf("在x=%f处的插值结果是: %f\n", xe, ye);  

    xe = 82.0;/* 填入要计算的x值 */;  
    ye = cubicSplineEval(x, y, b, c, d, N, xe);  
    printf("在x=%f处的插值结果是: %f\n", xe, ye);  

    return 0;  
}