finding coefficients for interpolation function

조회 수: 13 (최근 30일)
Melissa
Melissa 2013년 11월 25일
답변: Luke Von Neuman 2017년 6월 23일
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

답변 (3개)

Image Analyst
Image Analyst 2013년 11월 25일
"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
Image Analyst
Image Analyst 2013년 11월 25일
편집: Image Analyst 2013년 11월 25일
for k = 1 : length(coeff)
fprintf('coeff(%d) = %f\n', k, coeff(k));
end
Melissa
Melissa 2013년 11월 25일
Oh wow I didn't realize it was that simple. thank you!

댓글을 달려면 로그인하십시오.


Matt J
Matt J 2013년 11월 25일
편집: Matt J 2013년 11월 25일
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 2017년 6월 23일
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!

카테고리

Help CenterFile Exchange에서 Interpolation에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by