# finding coefficients for interpolation function

49 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

##### 0 Comments

### Answers (3)

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.

##### 0 Comments

Luke Von Neuman
on 23 Jun 2017

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!