Cubic Spline Interpolation Given End Conditions.

I'm writing a MATLAB program which accepts 3 inputs x (a vector containing the x values for interpolation), y (a vector containing the y values for interpolation) and a string specifying the type of cubic spline required ('natural', 'parabolically_terminated', 'not_a_knot') and then interpolates these points accordingly.

I've written a sample program (this clearly isn't the finished function) which computes the "a" vector (a vector containing the a coefficients that define the spline) of a natural cubic spline. The end conditions for a natural cubic spline are a1 = 0 and an = 0. The code is here:

x = [-2 -1 0 1 2];
y = [4 -1 2 1 8];
n=length(x) - 1; 
q = length(y) - 2;
h = x(2:n+1)-x(1:n);

v1=[h(1:n-1) 0];
v2=[1 2*(h(1:n-1)+h(2:n)) 1]; 
v3=[0 h(2:n)] ; 
M=diag(v1,-1)+diag(v2)+diag(v3,1);
d = [0 ((y(3: q + 2) - y(2: q + 1))./h(2:q + 1) - (y(2:q + 1) - y(1:q ))./h(1:q)) 0];
a = M\d'

I'm just wondering how I would alter this code to deal with a parabolically terminated cubic spline (which has end conditions -a_1 + a_2 = 0 and -a_{n - 1} + a_n = 0)? I know that the only changes will be to the first and last rows of M and d. I'm just not sure HOW they will differ from the natural cubic spline form. Thanks.

Answer 1

[Matt J]

I'm just wondering how I would alter this code to deal with a parabolically terminated cubic spline (which has end conditions -a_1 + a_2 = 0 and -a_{n - 1} + a_n = 0)?

Well, these are simply additional linear equalities, right? It looks like a matter of just adding them as additional rows to M and d'.

[David]

Actually, does the following code look correct? I've altered the first and last rows of A to allow for the end conditions, but I'm not so sure is it working properly.

n = length(x) - 1;                       %Prevents out of bounds errors.
h = x(2:n+1) - x(1:n);                   %Defines h.
if (strcmp(spline_type, 'natural') == 1) %If statements check which spline was chosen. 
    v1 = [h(1:n-1) 0];                   %Sub diag.
    v2 = [1 2 * (h(1:n-1)+h(2:n)) 1];    %Main diag.
    v3 = [0 h(2:n)] ;                    %Super diag.
    M = diag(v1,-1) + diag(v2) + diag(v3,1); %Creates a tridiagonal n x n sixe matrix.
elseif (strcmp(spline_type, 'parabolically_clamped') == 1)
    v1=[h(1:n-1) -1];                        %Notice differences in first and last rows to meet end conditions.
    v2=[-1 2*(h(1:n-1)+h(2:n)) 1]; 
    v3=[1 h(2:n)] ; 
    M=diag(v1,-1) + diag(v2) + diag(v3,1);   %Extra entries required to meet end conditions.
else
    v1=[h(1:n-1) 2];
    v2=[-1 2*h(1:n-1)+h(2:n) -1]; 
    v3=[2 h(2:n)] ;
    M = diag(v1,-1) + diag(v2) + diag(v3,1);
    M(1,3) = -1;
    M(n+1, n-1) = -1;                        %Extra entries required to meet end conditions.
end
d = [0 ((y(3: n + 1) - y(2: n ))./h(2:n) - (y(2:n) - y(1:n - 1 ))./h(1:n - 1)) 0]; %Defined the righ hand side matrix.
a = M\d';

[Matt J]

I can't tell. What I was picturing was that you would start with the M,d from your original post (assuming that's correct) and then simply add 2 additional rows at the end

[m,n]=size(M);

M(m+2,end-1:end)=[-1,1];
M(m+1,[1,2])=[-1,1];

d(end+2)=0;

Answer 2

[Matt J]

One other thing to consider is to use SPARSE to construct your M, since apparently, it is mostly tridiagonal.

This FEX submission (See in particular Example1.m) might offer a quick starting point:

http://www.mathworks.com/matlabcentral/fileexchange/26292-regular-control-point-interpolation-matrix-with-boundary-conditions