Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » Software » comp.soft-sys.matlab

Topic: Cubic Spline Interpolation
Replies: 1   Last Post: May 9, 2013 5:44 AM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Ben Vijayanayagam

Posts: 2
Registered: 5/9/13
Cubic Spline Interpolation
Posted: May 9, 2013 5:09 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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-')



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.