2011-01-20 8 views
2

私は正しく動作するようにMATLABコードを取得するのに問題があります!私は私に補間された多項式を与えるために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); 

クランプモードのための完全なコード:それは私だけ最初の2列を与え、このいずれかの

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つの仕事をすることができます知っている?

+1

複数の質問が含まれているようなコードを含め、膨大な量のコードがあります。一度に1つの質問に固執し、発生したエラーを再現するために必要な最小限の情報を使用してください。 –

答えて

-1

あなたのデータには3ポイントしかありません。キュービックに合うように4(またはそれ以上)のポイントが必要なので、おそらくアウトオブバウンドのエラーは、配列内の別のポイントを探しているコードから来ます。