0
#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;
}

Here's my part of source code written in visual studio 2017 environment. Globally initialized array x,y is given data. I have to do cubic spline interpolation with these data. But after execute, the spline results comes little differently at given points. I think there's no logical problem in spline and TDMA function. Maybe these problems are from storing irrational numbers. Can you guys explain why this error occur and how to fix the code?

hyojoon
  • 493
  • 1
  • 6
  • 17
  • 1
    https://stackoverflow.com/questions/588004/is-floating-point-math-broken – yano Apr 23 '18 at 15:48
  • Example: Assume `A[i] = (x[i] - x[i - 1])/6;` did `1.0/6`. Do you expect `A[i]` to have a value _exactly_ 1/6 or that `double A[i]` will take on the closest representable `double` to 1/6? `double` cannot represent every possible number - rounding occurs. The accumulation of the roundings on the various computations resulted in "results comes little differently at given points" that expected. For a general explanation see the above comment. Else append a specific question including your results and expectations. – chux - Reinstate Monica Apr 23 '18 at 18:41
  • Thanks to @chux, I completely understand – hyojoon Apr 24 '18 at 19:30

1 Answers1

0

I think your precision problem comes from using TDMA algorithm for the cubic spline interpolation. When doing cubic spline interpolation, you will need to solve a set of linear equations. TDMA algorithm is not the best way (precision wise) for solving linear equations. You can achieve a better accuracy by using alternatives such as GEPP (Gaussian Elimination with Partial Pivoting) or LU decomposition.

BTW, in your codes you should move the call to TDMA() out of cubicspline() function as you should only need to make the call once.

fang
  • 3,473
  • 1
  • 13
  • 19