c - 如何精确计算double类型数据来解决C中的三次样条插值？

Question

#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define SIZE 11

double cubicspline(double val);
void TDMA(void);

double x[SIZE] = { -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 };
double y[SIZE] = { 0.038, 0.058, 0.1, 0.2, 0.5, 1.0, 0.5, 0.2, 0.1, 0.058, 0.038 };
double M[SIZE];

void main(void) {
    FILE *fp1, *fp2;
    double val;
    int i;
    fp2 = fopen("cubicSpline_output.txt", "w");
    val = -1;
    while(val<=1.0) 
    {
        printf("------------------------------\n");
        printf("val : %lf\n",val);
        printf("spline result : %.12lf\n", cubicspline(val));
        fprintf(fp2, "%lf\t%lf\n", val, cubicspline(val));
        printf("------------------------------\n");
        val += 0.1;
    }
    fclose(fp2);
    system("pause");
}

void TDMA() {
    double A[SIZE-1], B[SIZE-1], C[SIZE-1], R[SIZE-1];
    double f[SIZE], g[SIZE];

    int i;

    for (i = 1; i < SIZE - 1; i++) 
    {
        A[i] = (x[i] - x[i - 1])/6;
        B[i] = (x[i + 1] - x[i - 1])/3;
        C[i] = (x[i + 1] - x[i])/6;
        R[i] = (y[i + 1] - y[i]) / 6 * C[i] - (y[i] - y[i - 1]) / 6 * A[i];

    }

    for (i = 1; i < SIZE - 2; i++) 
    {
        f[i + 1] = B[i + 1] - (C[i] * A[i + 1]) / B[i];
        g[i + 1] = R[i + 1] - (R[i] * A[i + 1]) / B[i];
    }

    M[0] = 0;
    M[SIZE - 1] = 0;
    M[SIZE - 2] = g[SIZE - 2] / f[SIZE - 2];

    for (i = SIZE - 3; i > 1; i--) 
    {
        M[i] = (g[i + 1] / f[i + 1]) - (C[i + 1]*M[i + 1]);
    }  

    M[1] = (R[1] - C[1]*M[2])/B[1];
}

double cubicspline(double val)
{
    int i, j;
    double result=0;

    TDMA();

    for (i = 1; i < SIZE; i++)
    {
        if (x[i] >= val) {
            j = i;
            break;
        }
    }
    result = ((pow((x[j] - val), 3.0) * M[j - 1]) + (pow((val - x[j - 1]), 3.0)) * M[j])/ 6 * (x[j] - x[j - 1]) + ((x[j] - val)*y[j - 1] + (val - x[j - 1])*y[j]) / (x[j] - x[j - 1]) - ((x[j] - x[j - 1])*((x[j] - val)*M[j - 1] + (val - x[j - 1])*M[j])) / 6;
    return result;
}

这是我在 Visual Studio 2017 环境下编写的部分源代码。全局初始化的数组 x,y 是给定的数据。我必须对这些数据进行三次样条插值。但执行后，样条结果在给定点处略有不同。我认为样条和TDMA函数没有逻辑问题。也许这些问题是由于存储无理数造成的。你们能解释一下为什么会出现这个错误以及如何修复代码吗？

Answer 1

我认为您的精度问题来自于使用 TDMA 算法进行三次样条插值。进行三次样条插值时，您需要求解一组线性方程。 TDMA 算法不是求解线性方程的最佳方法(精度方面)。您可以使用 GEPP(部分枢轴高斯消元法)或 LU 分解等替代方法来获得更高的精度。

顺便说一句，在您的代码中，您应该将对 TDMA() 的调用移出cubicspline() 函数，因为您只需要调用一次。