matlab中的三次樣條

Question

我無法讓matlab代碼正常工作！我在matlab中找到了一個三次樣條代碼給我插值多項式。和我簡單地給它的一個例子的工作：

僖= [0 0.05 0.1]

FI = [1 1.105171 1.221403]

連接」 = [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 

這一個只給我的第一個2列！

所以任何人都知道我可以如何讓這兩個工作？

Answer 1

您的數據只有3分。您需要4個（或更多）點來擬合立方，所以您的出界錯誤可能來自尋找數組中其他點的代碼。