MATLAB Answers

finding coefficients for interpolation function

34 views (last 30 days)
Melissa
Melissa on 25 Nov 2013
Answered: Luke Von Neuman on 23 Jun 2017
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

Answers (3)

Image Analyst
Image Analyst on 25 Nov 2013
"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.
  3 Comments
Melissa
Melissa on 25 Nov 2013
Oh wow I didn't realize it was that simple. thank you!

Sign in to comment.


Matt J
Matt J on 25 Nov 2013
Edited: Matt J on 25 Nov 2013
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
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.

Luke Von Neuman
Luke Von Neuman on 23 Jun 2017
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!

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!