matlab - matlab中的三次样条

Question

我无法让 matlab 代码正常工作!我在 matlab 中找到了一个三次样条代码来给我插值多项式。我只是给它一个工作的例子:

习 = [0 0.05 0.1]

Fi = [1 1.105171 1.221403]

Fi' = [2 _ 2.442806]

但它给了我这个错误:

??? 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);

这里是非结条件的完整代码

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

同样对于钳位条件，它给了我这个错误:

??? 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 );

钳位模式的完整代码:

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

对于这个，它只给我前两列!

所以有人知道我如何使这两个工作正常吗？

Answer 1

您的数据只有 3 分。您需要 4 个(或更多)点来拟合一个立方体，因此您的越界错误可能来自代码寻找数组中的另一个点。