2
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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

$\endgroup$
  • $\begingroup$ @Ozera: You're welcome :) $\endgroup$ – barak manos Apr 14 '14 at 22:08
  • $\begingroup$ 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. $\endgroup$ – Ozera Apr 14 '14 at 22:44
  • $\begingroup$ @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)$. $\endgroup$ – barak manos Apr 15 '14 at 17:03
  • $\begingroup$ @Ozera: What is the down-vote for????? Some of the people on this website are absolutely un@#$%believable!!! $\endgroup$ – barak manos Aug 30 '14 at 17:22
  • $\begingroup$ Odd to comment on a question from April, but I don't know why someone down voted you. $\endgroup$ – Ozera Sep 2 '14 at 8:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.