# Natural Cubic Spline Confusion

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.

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$.

• @Ozera: You're welcome :) – barak manos Apr 14 '14 at 22:08
• So i'm trying to do it by hand. So i'm following this example in my book i.imgur.com/o1rxB7C.png I calculated my interpolating conditions mathurl.com/lkyhm83 but i'm not sure how to do the derivative part. In particular, the continuity conditions. – Ozera Apr 14 '14 at 22:44
• @Ozera: The continuity and differentiability come along "naturally" with the way that $f''$ is defined. For the first and last points, $f''(x_0)=f''(x_n)=0$. For all other points, $f''(x_i)=\frac{3(h_{i+1}-h_i)}{(w_{i+1}+w_i)}$. Nevertheless, keep in mind that your final goal is $f(x)$ and not $f''(x)$. – barak manos Apr 15 '14 at 17:03
• @Ozera: What is the down-vote for????? Some of the people on this website are absolutely un@#\$%believable!!! – barak manos Aug 30 '14 at 17:22
• Odd to comment on a question from April, but I don't know why someone down voted you. – Ozera Sep 2 '14 at 8:00