Plotting a smooth curve from points
显示 更早的评论
I'm trying to plot a smooth line for 6 points in such a way that slope at each is zero. I have successfully achieved it for middle points using 'pchip' but I also want it for extreme points i.e. x=0 and x=15 which I am unable to do while using 'pchip'.
Here is my code
clear
clc
x = [0;3;6;9;12;15];
y = [0;1.9190;-3.2287;3.5133;-2.6825;1];
xi = linspace(min(x), max(x), 150); % Evenly-Spaced Interpolation Vector
yi = interp1(x, y, xi, 'pchip');
figure
plot(x, y, 'b')
hold on
plot(xi, yi, '-r')
hold off
grid
xlabel('X')
ylabel('Y')
legend('Original Data', 'Interpolation', 'Location', 'NE')
3 个评论
Image Analyst
2021-1-23
How do you know that the slope is exactly zero at the knot/training points? And why do you need the slope to be zero there?
Osama Anwar
2021-1-23
The slopes are very close to 0. Vertical lines show the original x-values without the endpoints.
x = [0;3;6;9;12;15];
y = [0;1.9190;-3.2287;3.5133;-2.6825;1];
xi = linspace(min(x), max(x), 150);
yi = interp1(x, y, xi, 'pchip');
% compute slopes
s = gradient(yi, 150);
plot(xi,s)
yline(0)
arrayfun(@xline, x(2:end-1)
采纳的回答
Bruno Luong
2021-1-24
编辑：Bruno Luong
2021-1-24
Direct analytic method using piecewise cublic polynomial. The curve is first-order differentiable, but not second order differentiable (as with sublic spline or result with my first answer using BSFK)
x = [0;3;6;9;12;15];
y = [0;1.9190;-3.2287;3.5133;-2.6825;1];
[xx,is] = sort(x(:));
yy = y(is);
yy = yy(:);
dx = diff(xx);
dy = diff(yy);
y0 = yy(1:end-1);
n = numel(xx)-1;
coefs = [-2*dy./(dx.^3), 3*dy./(dx.^2), 0*dy, y0];
pp = struct('form', 'pp',...
'breaks', xx(:)',...
'coefs', coefs,...
'pieces', n, ...
'order', 4,...
'dim', 1);
figure
xi = linspace(min(x),max(x));
yi = ppval(pp,xi);
plot(x,y,'b-o',xi,yi,'r');
xlim([min(x),max(x)])
grid on
Propably this is how John generates his curve
x = [0 1 1.5 2 4];
y = x;
10 个评论
Osama Anwar
2021-1-24
Bruno Luong
2021-1-24
Phi=[0,0,0,0,0;0.284629676546570,-0.830830026003772,1.30972146789057,-1.68250706566236,1.91898594722899;0.546200349457202,-1.08815592122522,0.372785597771793,1.39787738911579,-3.22870741511956;0.763521118433367,-0.594351144437141,-1.20361562377557,0.521108558113203,3.51333709166613;0.918985947228994,0.309721467890570,-0.715370323453430,-1.83083002600377,-2.68250706566236;1,1,1,1,1]
x = [-3;0;3;6;9;12;15;18];
y = [2;Phi(:,1);-2];
% ...
Osama Anwar
2021-1-24
编辑：Osama Anwar
2021-1-24
Osama Anwar
2021-1-24
Osama Anwar
2021-1-24
编辑：Osama Anwar
2021-1-24
Bruno Luong
2021-1-24
You meant this:
Simply reverse in the plot
Phi=[0,0,0,0,0;0.284629676546570,-0.830830026003772,1.30972146789057,-1.68250706566236,1.91898594722899;0.546200349457202,-1.08815592122522,0.372785597771793,1.39787738911579,-3.22870741511956;0.763521118433367,-0.594351144437141,-1.20361562377557,0.521108558113203,3.51333709166613;0.918985947228994,0.309721467890570,-0.715370323453430,-1.83083002600377,-2.68250706566236;1,1,1,1,1]
x = [0;3;6;9;12;15];
y = [Phi(:,1)];
[xx,is] = sort(x(:));
yy = y(is);
yy = yy(:);
dx = diff(xx);
dy = diff(yy);
y0 = yy(1:end-1);
n = numel(xx)-1;
coefs = [-2*dy./(dx.^3), 3*dy./(dx.^2), 0*dy, y0];
pp = struct('form', 'pp',...
'breaks', xx(:)',...
'coefs', coefs,...
'pieces', n, ...
'order', 4,...
'dim', 1);
figure
xi = linspace(min(x),max(x),257);
yi = ppval(pp,xi);
plot(y,x,'bo',yi,xi,'r');
grid on
xlabel('y')
ylabel('x')
Osama Anwar
2021-1-24
Osama Anwar
2021-1-27
Bruno Luong
2021-1-27
Sorry you are on your own if you use symbolic variables. Once you modify/adap the code it's yours.
I don't have this toolbox to help, even if I want to .
Osama Anwar
2021-1-27
更多回答（2 个）
You could add values to the beginning and end to make the curve continuing in both directions. The example below uses hard-coded values but it wouldn't be difficult to programmatically extend the trend in the opposite y-direction on each side of X.
x = [-3;0;3;6;9;12;15;18];
% ^^ ^^ added
y = [2;0;1.9190;-3.2287;3.5133;-2.6825;1;-2];
% ^ ^^ added
xi = linspace(min(x), max(x), 150);
yi = interp1(x, y, xi, 'pchip');
% Trim the additions
rmIdx = xi<0 | xi>15;
xi(rmIdx) = [];
yi(rmIdx) = [];
x([1,end]) = [];
y([1,end]) = [];
figure
plot(x, y, 'b')
hold on
plot(xi, yi, '-r')
hold off
grid
xlabel('X')
ylabel('Y')
legend('Original Data', 'Interpolation', 'Location', 'NE')
13 个评论
Osama Anwar
2021-1-23
Adam Danz
2021-1-23
I just update the answer to add the 5 lines under "Trim the additions".
To make the solution programmatic you just need to implement these steps.
- Determine if the slope of + or - at the end points. That's easy. At each endpoint you just need to determine if the y-value of the end point is greater or less than the y-value of the neighboring coordinate.
- Compute the median x-interval: median(diff(x))
- Compute the median y-distance from center: median(abs(y)) (since your curve is centered at y=0)
- Add a coordinate to the right of the curve by adding the median x-interval to the x-value of the endpoint and set the y values as the median y-distance with its sign opposite to the end point's slope.
- Do the same for the left end-point by subtracting the median x-interval.
Osama Anwar
2021-1-23
Adam Danz
2021-1-23
That's true. x-interval isn't important.
Osama Anwar
2021-1-23
Adam Danz
2021-1-23
Yes, that's true, too. Thanks for following up!
John D'Errico
2021-1-23
I can think of a few ways to do this. Adam's solution is a nice one though. +1
Hmm, as I think and re-read the question, this solution works ONLY if each point is an extremeum.
Is the target to have an interpolant that has zero slope at each point, even for a sequence of points that follow a straight line? For example, consider the points:
x = [0 1 1.5 2 4];
y = x;
plot(x,y,'-o')
Or, is this next curve your goal?
Adam's solution will not produce the latter curve in general. So which is your goal?
Osama Anwar
2021-1-24
Osama Anwar
2021-1-24
Osama Anwar
2021-1-24
编辑：Osama Anwar
2021-1-24
Adam Danz
2021-1-24
Yeah, it's not a robust hack.
I'm following the conversation. John and Bruno have great points.
Osama Anwar
2021-1-24
Bruno Luong
2021-1-24
编辑：Bruno Luong
2021-1-24
No extra points needed (but you might add to twist the shape of the curve in the first and last interval),
Spline order >= 8th is needed using my FEX
x = [0;3;6;9;12;15];
y = [0;1.9190;-3.2287;3.5133;-2.6825;1];
interp = struct('p', 0, 'x', x, 'v', y);
slope0 = struct('p', 1, 'x', x, 'v', 0*x);
% Download BSFK function here
% https://www.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation
pp = BSFK(x,y, 9, length(x)-1, [], struct('KnotRemoval', 'none', 'pntcon', [interp slope0]));
% Check
figure
xi = linspace(min(x),max(x));
yi = ppval(pp,xi);
plot(x,y,'b',xi,yi,'r');
for xb=pp.breaks
xline(xb);
end
grid on
4 个评论
Bruno Luong
2021-1-24
With Adam's trick
x = [0;3;6;9;12;15];
y = [0;1.9190;-3.2287;3.5133;-2.6825;1];
% Adam's trick
xx = [-3;x;18];
yy = [2;y;-2];
interp = struct('p', 0, 'x', x', 'v', y');
slope0 = struct('p', 1, 'x', x, 'v', 0*x);
% Download BSFK function here
% https://www.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation
pp = BSFK(xx,yy, 9, length(x)-1, [], struct('KnotRemoval', 'none', 'pntcon', [interp slope0]));
% Check
figure
xi = linspace(min(x),max(x));
yi = ppval(pp,xi);
plot(x,y,'b',xi,yi,'r');
for xb=pp.breaks
xline(xb);
end
xlim([min(x),max(x)])
grid on
Bruno Luong
2021-1-24
Monotonic data
x = [0;3;6;9;12;15];
y = x;
% Adam's trick
xx = [-3;x;18];
yy = xx;
interp = struct('p', 0, 'x', x', 'v', y');
slope0 = struct('p', 1, 'x', x, 'v', 0*x);
% Download BSFK function here
% https://www.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation
pp = BSFK(xx,yy, 9, length(x)-1, [], struct('KnotRemoval', 'none', 'pntcon', [interp slope0]));
% Check
figure
xi = linspace(min(x),max(x));
yi = ppval(pp,xi);
plot(x,y,'b-o',xi,yi,'r');
for xb=pp.breaks
xline(xb);
end
xlim([min(x),max(x)])
grid on
Osama Anwar
2021-1-24
Bruno Luong
2021-1-25
编辑：Bruno Luong
2021-1-26
The shape might not meet your (un-written) expectation but it definitively meets every requiremenrs you state in the question. The interpolating solution is trully "smooth" (I believe up to the 7th derivative order is continue) with zero slope at the give abscissa.
This illustres one of the difficulty using spline fitting/interpolating: the chosen boundary conditions affects globally the shape of the interpolation.
It might be still useful for futur readers.
类别
在 帮助中心 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
