Natural Cubic Spline Confusion

Question

Find the natural cubic spline which interpolates the data points $(1,0),\; (2,1),\; (3,0), \; (4,1), \; (5,0) $.

I know how to check if a piecewise function is a natural cubic spline, but I don't really know how to find a function that interpolates data points like that.

Answer 1

Here is a cubic-spline interpolation for the $5$ points given in your question:

$ f(x)= \begin{cases} -0.5(x-1)^3 + 1.5(x-1) & \text{$ 1 \leq x \leq 2$}\\ (x-2)^3 - 1.5(x-2)^2 - 0.5(x-2) + 1 & \text{$ 2 \leq x \leq 3$}\\ -(x-3)^3 + 1.5(x-3)^2 + 0.5(x-3) & \text{$ 3 \leq x \leq 4$}\\ 0.5(x-4)^3 - 1.5(x-4)^2 + 1 & \text{$ 4 \leq x \leq 5$}\\ \end{cases} $

Here is a piece of C code for any given number of points $(x_0,y_0),(x_1,y_1),\ldots,(x_N,y_N)$:

void Spline(double x[N+1],double y[N+1], // input
            double A[N],double B[N],     // output
            double C[N],double D[N])     // output
{
    double w[N];
    double h[N];
    double ftt[N+1];

    for (int i=0; i<N; i++)
    {
        w[i] = (x[i+1]-x[i]);
        h[i] = (y[i+1]-y[i])/w[i];
    }

    ftt[0] = 0;
    for (int i=0; i<N-1; i++)
        ftt[i+1] = 3*(h[i+1]-h[i])/(w[i+1]+w[i]);
    ftt[N] = 0;

    for (int i=0; i<N; i++)
    {
        A[i] = (ftt[i+1]-ftt[i])/(6*w[i]);
        B[i] = ftt[i]/2;
        C[i] = h[i]-w[i]*(ftt[i+1]+2*ftt[i])/6;
        D[i] = y[i];
    }
}

void PrintSpline(double x[N+1],           // input
                 double A[N],double B[N], // input
                 double C[N],double D[N]) // input
{
    for (int i=0; i<N; i++)
    {
        printf("%f <= x <= %f : f(x) = ",x[i],x[i+1]);
        printf("%f(x-%f)^3 + ",A[i],x[i]);
        printf("%f(x-%f)^2 + ",B[i],x[i]);
        printf("%f(x-%f) + "  ,C[i],x[i]);
        printf("%f\n"         ,D[i]);
    }
}

Here is a piece of Python code for any given number of points $(x_0,y_0),(x_1,y_1),\ldots,(x_N,y_N)$:

class Point:
    def __init__(self,x,y):
        self.x = 1.0*x
        self.y = 1.0*y

def Spline(points):
    N   = len(points)-1
    w   =     [(points[i+1].x-points[i].x)      for i in range(0,N)]
    h   =     [(points[i+1].y-points[i].y)/w[i] for i in range(0,N)]
    ftt = [0]+[3*(h[i+1]-h[i])/(w[i+1]+w[i])    for i in range(0,N-1)]+[0]
    A   =     [(ftt[i+1]-ftt[i])/(6*w[i])       for i in range(0,N)]
    B   =     [ftt[i]/2                         for i in range(0,N)]
    C   =     [h[i]-w[i]*(ftt[i+1]+2*ftt[i])/6  for i in range(0,N)]
    D   =     [points[i].y                      for i in range(0,N)]
    return A,B,C,D

def PrintSpline(points,A,B,C,D):
    for i in range(0,len(points)-1):
        func = str(points[i].x)+' <= x <= '+str(points[i+1].x)+' : f(x) = '
        components = []
        if A[i]:
            components.append(str(A[i])+'(x-'+str(points[i].x)+')^3')
        if B[i]:
            components.append(str(B[i])+'(x-'+str(points[i].x)+')^2')
        if C[i]:
            components.append(str(C[i])+'(x-'+str(points[i].x)+')')
        if D[i]:
            components.append(str(D[i]))
        if components:
            func += components[0]
            for i in range (1,len(components)):
                if components[i][0] == '-':
                    func += ' - '+components[i][1:]
                else:
                    func += ' + '+components[i]
            print func
        else:
            print func+'0'

def Example():
    points = [Point(1,0),Point(2,1),Point(3,0),Point(4,1),Point(5,0)]
    A,B,C,D = Spline(points)
    PrintSpline(points,A,B,C,D)

Please note that the two pieces of code above assume $x_0 < x_1 < \ldots < x_N$.