34 views (last 30 days)

Show older comments

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

Image Analyst
on 25 Nov 2013

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
on 23 Jun 2017

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

Start Hunting!