i have trouble getting a matlab code to work properly! i found a cubic spline code in matlab to give me the interpolated polynomial. and i simply give it an example to work:
Xi = [0 0.05 0.1]
Fi = [1 1.105171 1.221403]
Fi' = [2 _ 2.442806]
but it gives me this error:
??? Attempted to access du(1); index out of bounds because numel(du)=0.
Error in ==> cubic_nak at 53
du(1) = du(1) - hi(1)^2 / hi(2);
here is the full code for not a knot condition
function csn = cubic_nak ( xi, fi )
%CUBIC_NAK compute the cubic spline interpolant, subject to
% "not-a-knot" boundary conditions, associated with a
% given set of interpolating points and function values
%
% calling sequences:
% csn = cubic_nak ( xi, fi )
% cubic_nak ( xi, fi )
%
% inputs:
% xi vector containing the interpolating points
% (must be sorted in ascending order)
% fi vector containing function values
% the i-th entry in this vector is the function
% value associated with the i-th entry in the 'xi'
% vector
%
% output:
% csn five column matrix containing the information which
% defines the "not-a-knot" cubic spline interpolant
% - first column: interpolating points
% - second column: function values
% - third column: coefficients of linear terms
% - fourth column: coefficients of quadratic terms
% - fifth column: coefficients of cubic terms
%
% NOTE:
% to evaluate the "not-a-knot" cubic spline interpolant apply
% the routine SPLINE_EVAL to the output of this routine
%
n = length ( xi );
m = length ( fi );
if ( n ~= m )
disp ( 'number of ordinates and number of function values must be equal' )
return
end
for i = 1 : n-1
hi(i) = xi(i+1) - xi(i);
end
for i = 1 : n-2
dd(i) = 2.0 * ( hi(i) + hi(i+1) );
ri(i) = (3.0/hi(i+1))*(fi(i+2)-fi(i+1))-(3.0/hi(i))*(fi(i+1)-fi(i));
end
dd(1) = dd(1) + hi(1) + hi(1)^2 / hi(2);
dd(n-2) = dd(n-2) + hi(n-1) + hi(n-1)^2 / hi(n-2);
du = hi(2:n-2);
dl = du;
du(1) = du(1) - hi(1)^2 / hi(2);
dl(n-3) = dl(n-3) - hi(n-1)^2 / hi(n-2);
temp = tridiagonal ( dl, dd, du, ri );
c = zeros ( n,1 );
d = c; b = c;
c(2:n-1) = temp;
c(1) = ( 1 + hi(1) / hi(2) ) * c(2) - hi(1) / hi(2) * c(3);
c(n) = ( 1 + hi(n-1) / hi(n-2) ) * c(n-1) - hi(n-1) / hi(n-2) * c(n-2);
for i = 1 : n-1
d(i) = (c(i+1)-c(i))/(3.0*hi(i));
b(i) = (fi(i+1)-fi(i))/hi(i) - hi(i)*(c(i+1)+2.0*c(i))/3.0;
end
if ( nargout == 0 )
disp ( [ xi' fi' b c d ] )
else
csn = [ xi' fi' b c d ];
end
also for the clamped condition it gives me this error:
??? Undefined function or method 'tridiagonal' for input arguments of type 'double'.
Error in ==> cubic_clamped at 55
c = tridiagonal ( hi(1:n-1), dd, hi(1:n-1), ri );
??? Input argument "xi" is undefined.
Error in ==> cubic_clamped at 35
n = length ( xi );
the full code for clamped mode:
function csc = cubic_clamped ( xi, fi, fpa, fpb )
%CUBIC_CLAMPED compute the cubic spline interpolant, subject to
% "clamped" boundary conditions, associated with a
% given set of interpolating points and function values
%
% calling sequences:
% csc = cubic_clamped ( xi, fi, fpa, fpb )
% cubic_clamped ( xi, fi, fpa, fpb )
%
% inputs:
% xi vector containing the interpolating points
% (must be sorted in ascending order)
% fi vector containing function values
% the i-th entry in this vector is the function
% value associated with the i-th entry in the 'xi'
% vector
% fpa derivative value at left endpoint; i.e., xi(1)
% fpb derivative value at right endpoint; i.e., xi(n)
%
% output:
% csn five column matrix containing the information which
% defines the "clamped" cubic spline interpolant
% - first column: interpolating points
% - second column: function values
% - third column: coefficients of linear terms
% - fourth column: coefficients of quadratic terms
% - fifth column: coefficients of cubic terms
%
% NOTE:
% to evaluate the "clamped" cubic spline interpolant apply
% the routine SPLINE_EVAL to the output of this routine
%
n = length ( xi );
m = length ( fi );
if ( n ~= m )
disp ( 'number of ordinates and number of function values must be equal' )
return
end
for i = 1 : n-1
hi(i) = xi(i+1) - xi(i);
end
dd(1) = 2.0*hi(1); dd(n) = 2.0*hi(n-1);
ri(1) = (3.0/hi(1))*(fi(2)-fi(1)) - 3.0 * fpa;
ri(n) = 3.0 * fpb - (3.0/hi(n-1))*(fi(n)-fi(n-1));
for i = 1 : n-2
dd(i+1) = 2.0 * ( hi(i) + hi(i+1) );
ri(i+1) = (3.0/hi(i+1))*(fi(i+2)-fi(i+1))-(3.0/hi(i))*(fi(i+1)-fi(i));
end
disp ( [dd' ri'] )
c = tridiagonal ( hi(1:n-1), dd, hi(1:n-1), ri );
d = zeros ( n,1 );
b = d;
for i = 1 : n-1
d(i) = (c(i+1)-c(i))/(3.0*hi(i));
b(i) = (fi(i+1)-fi(i))/hi(i) - hi(i)*(c(i+1)+2.0*c(i))/3.0;
end
if ( nargout == 0 )
disp ( [ xi' fi' b c' d ] )
else
csc = [ xi' fi' b c' d ];
end
for this one it only gives me the first 2 columns!
so anyone knows how i can make these 2 work?
Answer
Your data only has 3 points. You need 4 (or more) points to fit a cubic, so your out-of-bounds error probably comes from the code looking for another point in the array.