Natural Cubic Spline Interpolation in C

spline.c

/** Numerical Analysis 9th ed - Burden, Faires (Ch. 3 Natural Cubic Spline, Pg. 149) */
#include <stdio.h>

int main() {
    /** Step 0 */
    int n, i, j;
    scanf("%d", &n);
    n--;
    float x[n + 1], a[n + 1], h[n], A[n], l[n + 1],
        u[n + 1], z[n + 1], c[n + 1], b[n], d[n];
    for (i = 0; i < n + 1; ++i) scanf("%f", &x[i]);
    for (i = 0; i < n + 1; ++i) scanf("%f", &a[i]);

    /** Step 1 */
    for (i = 0; i <= n - 1; ++i) h[i] = x[i + 1] - x[i];

    /** Step 2 */
    for (i = 1; i <= n - 1; ++i)
        A[i] = 3 * (a[i + 1] - a[i]) / h[i] - 3 * (a[i] - a[i - 1]) / h[i - 1];

    /** Step 3 */
    l[0] = 1;
    u[0] = 0;
    z[0] = 0;

    /** Step 4 */
    for (i = 1; i <= n - 1; ++i) {
        l[i] = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * u[i - 1];
        u[i] = h[i] / l[i];
        z[i] = (A[i] - h[i - 1] * z[i - 1]) / l[i];
    }

    /** Step 5 */
    l[n] = 1;
    z[n] = 0;
    c[n] = 0;

    /** Step 6 */
    for (j = n - 1; j >= 0; --j) {
        c[j] = z[j] - u[j] * c[j + 1];
        b[j] = (a[j + 1] - a[j]) / h[j] - h[j] * (c[j + 1] + 2 * c[j]) / 3;
        d[j] = (c[j + 1] - c[j]) / (3 * h[j]);
    }

    /** Step 7 */
    printf("%2s %8s %8s %8s %8s\n", "i", "ai", "bi", "ci", "di");
    for (i = 0; i < n; ++i)
        printf("%2d %8.2f %8.2f %8.2f %8.2f\n", i, a[i], b[i], c[i], d[i]);
    return 0;
}

Link to the book

http://ins.sjtu.edu.cn/people/mtang/textbook.pdf#page=167
https://fac.ksu.edu.sa/sites/default/files/numerical_analysis_9th.pdf#page=167

Input Format

n
x0 x1 x2 ... xn
y0 y1 y2 ... yn

where
n = no. of points
xi = abscissa of point i
yi = ordinate of point i

Output Format

i ai bi ci di
...
...

where
i = ith spline polynomial
a, b, c, d = coefficients of the polynomial

Thank you, thank you so much, so very useful and I learned how the spline works using its coefficients, thank you again 👍
NOTE!!!
You have an error in Line 28, you wrote:
l[i] = 2 * (x[i + 1] - x[i]) - h[i - 1] * u[i - 1];

you forgot the minus 1 (-1), the correct form is this:
l[i] = 2 * (x[i + 1] - x[i-1]) - h[i - 1] * u[i - 1];

Link to the book is broken. I found this one working: http://fac.ksu.edu.sa/sites/default/files/textbook-9th_edition.pdf

@ekalkan
thank you for the link.

@ivanfe639
👍

Thanks @ivanfe639 and @ekalkan, corrected the code in 28th line and updated the book link.

Thanks a lot for the code,
I modified it so that it can interpolate an array of 'n' values to an array of 'x' values.

#include <stdio.h>
#define CURRENT_NUMBER_OF_POINTS 40
#define NUMBER_INTERPOLATED_POINTS 2000
int main() {
/** Step 0 */
int n=CURRENT_NUMBER_OF_POINTS, i, j;

n--;
float x[n + 1], a[n + 1], h[n], A[n], l[n + 1],
        u[n + 1], z[n + 1], c[n + 1], b[n], d[n];
printf("enter values of points\n");
for (i = 0; i < n + 1; ++i) scanf("%f", &a[i]);


for (i = 0; i < n + 1; ++i) x[i]=i;

/** Step 1 */
for (i = 0; i <= n - 1; ++i) h[i] = x[i + 1] - x[i];

/** Step 2 */
for (i = 1; i <= n - 1; ++i)
    A[i] = 3 * (a[i + 1] - a[i]) / h[i] - 3 * (a[i] - a[i - 1]) / h[i - 1];

/** Step 3 */
l[0] = 1;
u[0] = 0;
z[0] = 0;

/** Step 4 */
for (i = 1; i <= n - 1; ++i) {
    l[i] = 2 * (x[i + 1] - x[i-1]) - h[i - 1] * u[i - 1];
    u[i] = h[i] / l[i];
    z[i] = (A[i] - h[i - 1] * z[i - 1]) / l[i];
}

/** Step 5 */
l[n] = 1;
z[n] = 0;
c[n] = 0;

/** Step 6 */
for (j = n - 1; j >= 0; --j) {
    c[j] = z[j] - u[j] * c[j + 1];
    b[j] = (a[j + 1] - a[j]) / h[j] - h[j] * (c[j + 1] + 2 * c[j]) / 3;
    d[j] = (c[j + 1] - c[j]) / (3 * h[j]);
}

/** Step 7 */
float interpolated_x[NUMBER_INTERPOLATED_POINTS];

/* Evaluate cubic spline equation at interpolated x values*/
for (i = 0; i < NUMBER_INTERPOLATED_POINTS; ++i) {
float x_val = (float)i * n / (float)(NUMBER_INTERPOLATED_POINTS - 1);
int j = (int)x_val;
float dx = x_val - j;
interpolated_x[i] = a[j] + b[j]dx + c[j]dxdx + d[j]dxdxdx;
}

    for (int i = 0; i < NUMBER_INTERPOLATED_POINTS; i++) {
        printf( "%d == %f \n", i,interpolated_x[i]);
    }
 
return 0;

}

Hey so im making a hyper specific library for matlab, to use in a project, is it fine if i use this? and if i do then how should i credit you?

after testing it a bunch, if you want to use the values in a polynomial to find out querry points, a , b , c and d are in the reverse order,

i compared it to here : https://www.bragitoff.com/2018/02/cubic-spline-piecewise-interpolation-c-program/?unapproved=60552&moderation-hash=d514da745560ff352ee385453b9346e1#comment-60552

Hey so im making a hyper specific library for matlab, to use in a project, is it fine if i use this? and if i do then how should i credit you?

@Pin09091 yes of course. Feel free to use it, mention link to the gist of my name. I'm glad this code is useful for you