finding coefficients for interpolation function

Good Morning All,

I have created a function called interpolation which gives the user an option to choose between linear, cubic spline and akima interpolation techniques when chosen. Currently the interpolation values are given and I would like to have the coefficients of the function displayed as well. Does anyone have any suggestions on how to do this?

I am trying to obtain the coefficients so I can then take the derivative of said polynomial and use this for a newton raphson method. So the function and its derivative are vital to complete my task.

Thanks for all of your help,

Mel

function yi = interpolation(x,y,xi,method)

%function yi = interpolation(x,y,xi,method) is an interpolation function
%designed for force-displacement curve. The user has the option to select
%between linear, cubic spline and akima for methods of evaluation.

% User Input
%   x: interpolation nodes [nx1]
%   y: interpolation values [nx1]
%   xi: evaluation points for the inerpolant
%   method: specifies alternate methods for interpolation
%           'linear'   - linear interpolation
%           'cspline'   - piecewise cubic spline interpolation 
%           'akima'   - akima spline interpolation 
% Output
%   yi: desired interpolation values [nx1]

%Switch cases for method:
switch method

%Linear Case    
    case 'linear'
n = length(x)-1;
%Initialize Output
yi = NaN*zeros(size(xi));
% Piecewise Slopes
m = diff(y) ./ diff(x);
%Solving for yi
for k=1:n
  %Find points within xi interval to evaluate  
  xii = find( xi>=x(k) & xi<=x(k+1) ); 
  yi(xii) = y(k) + m(k)*(xi(xii)-x(k));  
end
yi=yi';

%Cubic Case
    case 'cspline'
%Ensure column vectors       
x = x(:);  y = y(:);
n = length(x)-1;
m = diff(x);  
ddx = diag(m);  ddx2 = diag(m.^2);  ddx3 = diag(m.^3);
D = toeplitz( [1;zeros(n-2,1)], [1 -1 zeros(1,n-2)] );

% Interpolation conditions (n rows)
A1 = [ ddx3 ddx2 ddx ];
rhs1 = diff(y);

% Continuity of yi' and yi'' (n-1 rows each)
A2 = [ 3*ddx2(1:n-1,:) 2*ddx(1:n-1,:) D ];
rhs2 = zeros(n-1,1);
A3 = [ 3*ddx(1:n-1,:) D zeros(n-1,n) ];
rhs3 = zeros(n-1,1);

% Not-a-knot end conditions (2 rows)
NAK = [ [1 -1 zeros(1,3*n-2)]; [zeros(1,n-2) 1 -1 zeros(1,2*n)] ];
rhs4 = [0;0];

% Assemble and solve system
coeff = [ A1;A2;A3;NAK ] \ [rhs1;rhs2;rhs3;rhs4];
coeff = reshape( coeff, [n,3] );
coeff(:,4) = y(1:n);  % known constant term in cubic

yi = zeros(size(xi));
%Solving for yi
for k=1:n
%Find points within xi interval to evaluate  
  xii = (xi>=x(k)) & (xi<=x(k+1));
  yi(xii) = polyval( coeff(k,:), xi(xii)-x(k) )';
end
yi=yi';

%Case Akima
    case'akima'
%Ensure column vectors  
x=x(:); y=y(:); xi=xi(:); n=length(x);
dx=diff(x);
m=diff(y)./dx;
mm=2*m(1)-m(2);     mmm=2*mm-m(1);     % augment at left
mp=2*m(n-1)-m(n-2); mpp=2*mp-m(n-1);   % augment at right
m1=[mmm; mm; m; mp; mpp];              % slopes
dm=abs(diff(m1)); f1=dm(3:n+2); f2=dm(1:n); f12=f1+f2;
id=find(f12 > 1e-8*max(f12)); b=m1(2:n+1);
b(id)=(f1(id).*m1(id+1)+f2(id).*m1(id+2))./f12(id);
c=(3*m-2*b(1:n-1)-b(2:n))./m;
d=(b(1:n-1)+b(2:n)-2*m)./m.^2;

[ncnt,bin]=histc(xi,x);
bin=min(bin,n-1);
bb=bin(1:length(xi));
wj=xi-x(bb);
yi=((wj.*d(bb) +c(bb)).*wj+b(bb)).*wj+y(bb);
end

Answer 1

[Image Analyst]

"I would like to have the coefficients of the function displayed" - I would either just put the name of the variable all by itself on a line of code without a semicolon at the end of the line, OR use fprintf() to get more precise control of exactly what is displayed.

[Melissa]

Thanks ImageAnalyst for the response. My concern is its solving for every iteration so I have those results for every "k" how would I list the actually coefficients or use fprint?

[Image Analyst]

for k = 1 : length(coeff)
  fprintf('coeff(%d) = %f\n', k, coeff(k));
end

[Melissa]

Oh wow I didn't realize it was that simple. thank you!

Answer 2

[Matt J]

I know nothing about Akima interpolation. Linear interpolation uses non-differentiable piecewise linear polynomials, so Newton Raphson won't be applicable with that. Since you are dealing with 1D interpolation, I wonder why fzero() wouldn't be better for whatever you are trying to do.

Regardless, here is a FEX tool for taking the derivatives of differentiable splines

http://www.mathworks.com/matlabcentral/fileexchange/24996-spline-derivative

You can also use the interp1 syntax

 pp = interp1(x,Y,method,'pp')

which returns the piece-wise polynomial description of the interpolation. The pp object contains coefficient information which can be extracted with UNMKPP.

Answer 3

[Luke Von Neuman]

Hello! Hello guys! I've tried to understand the part of code, which describes 'cspline' method and transform it in analytic model. I can't understand exactly how it's implemented. First of all I believe that it's linked to this algorithm from wiki, ---> https://en.wikipedia.org/wiki/Spline_interpolation, but doesn't match.Could you help me to transform it in mathematical model? Al least, a link which teach me better. I need this urgently. Thanks a lot!