finding coefficients for interpolation function
8 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(3 个)
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
2013-11-25
编辑:Image Analyst
2013-11-25
for k = 1 : length(coeff)
fprintf('coeff(%d) = %f\n', k, coeff(k));
end
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.
0 个评论
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!
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!