Thread Subject:

manually creating a spline

Hi, I am trying to manually create a natural cubic spline in Matlab without using the built in spline function. The test problem has 3 spans (4 control points), and I am keeping it general enough so that it can handle unequally spaced points.

Ultimately, I want to extract interpolated values within each span but am having difficulty in getting the code right. Can anyone with expertise in this area shed some light on what I'm doing wrong with the below code please? Thanks a lot.


%program to simulate spline behaviour for 3 spans
% Control point coordinates
x(1) = 1;
x(2) = 6;
x(3) = 11;
x(4) = 16;

y(1) = 1;
y(2) = 3;
y(3) = 2.5;
y(4) = 0.5;

%spans
h1 = x(2) - x(1);
h2 = x(3) - x(2);
h3 = x(4) - x(3);

%global matrix for second derivatives
B = 6*[0;
    ((y(3)-y(2))/h2)-((y(2)-y(1))/h1);
    ((y(4)-y(3))/h3)-((y(3)-y(2))/h2);
    0]

A = [1 0 0 0;
    h1 2*(h1+h2) h2 0;
    0 h2 2*(h2+h3) h3;
    0 0 0 1]

% constants
reqs = inv(A)*B
cst1 = reqs(2,1)
cst2 = reqs(3,1)

% coefficients
a1 = (cst1-0)/(6*h1);
a2 = (cst2-cst1)/(6*h2);
a3 = (0-cst2)/(6*h3);

b1 = 0;
b2 = cst1/2;
b3 = cst2/2;

c1 = (y(2)-y(1))/h1 - cst1*h1/6 - 0;
c2 = (y(3)-y(2))/h2 - cst2*h2/6 - cst1*h2/3;
c3 = (y(4)-y(3))/h3 - 0 - cst2*h3/3;

d1 = y(1);
d2 = y(2);
d3 = y(3);
d4 = y(4);

% spline coeffs
A1=[a1; b1; c1; d1];
A2=[a2; b2; c2; d2];
A3=[a3; b3; c3; d3];

% interpolated positions
t1=1:6;
t2=6:11;
t3=11:16;

% splines
spl1=(a1*t1.^3)+(b1*t1.^2)+(c1*t1)+d1
spl2=(a2*t2.^3)+(b2*t2.^2)+(c2*t2)+d2
spl3=(a3*t3.^3)+(b3*t3.^2)+(c3*t3)+d3

plot(t1,spl1,'b-',t2,spl2,'b-',t3,spl3,'b-')