Fitting a monotonically increasing spline function

2022 9 5
1 个回答
回答已采纳
更新时间：2022 9 6
25 次查看（30 天）
Ran in:

0 个投票

Ran in:
I want to fit a monotonously increasing smooth spline function for a dataset
x = [0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 20.85, 21.]
x = 1×17
0 0.7500 1.8000 2.2500 3.7500 4.5000 6.4500 6.7500 7.5000 8.3250 10.8750 11.2500 12.5250 12.7500 15.0000 20.8500 21.0000
y = [2.83811035, 2.81541896, 3.14311655, 3.22373554, 3.43033456, 3.50433385, 3.66794514, 3.462296, 3.59480959,
3.56250726, 3.6209845, 3.63034523, 3.68238915, 3.69096892, 3.75560395, 3.83545191, 3.90419498]
Error using vertcat
Dimensions of arrays being concatenated are not consistent.
The current fit using interp1d looks like the above. I would like to know how to fit a monotonously increasing spline function.

请先登录，再回答此问题。

 采纳的回答

Bruno Luong
Bruno Luong 2022-9-5

2 个投票

One way is to use my https://uk.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation
x = [0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 20.85, 21.]
y = [2.83811035, 2.81541896, 3.14311655, 3.22373554, 3.43033456, 3.50433385, 3.66794514, 3.462296, 3.59480959,3.56250726, 3.6209845, 3.63034523, 3.68238915, 3.69096892, 3.75560395, 3.83545191, 3.90419498]
nknots=10;
opt=struct('shape',struct('p',1,'lo',zeros(1,nknots),'up',inf(1,nknots)));
pp=BSFK(x,y,4,nknots,[],opt); %FEX
xi=linspace(min(x),max(x),100);
yi=ppval(pp,xi);
plot(xi,yi,'-',x,y,'or')

7 个评论
显示 5更早的评论 隐藏 5更早的评论

Deepa Maheshvare
Deepa Maheshvare 2022-9-6
@Bruno Luong Thanks so much. This is great.
Could you please explain a bit about how nknots=10 is set?
Also, may I know if there is a python implementation of this function.
Bruno Luong
Bruno Luong 2022-9-6
I select 10 knots because you have 13 points, do 10 should be enough. I believe at tte end BSFK remove the knots and keep relevent 4 or 5.
No there is no equivalent python implementation of my function.
Deepa Maheshvare
Deepa Maheshvare 2022-9-6
@Bruno Luong, Could we get the polynomial function of the fitted curve?
Bruno Luong
Bruno Luong 2022-9-6
You could use the same technique I gave here :
https://uk.mathworks.com/matlabcentral/answers/1794985-how-to-constrain-the-resulting-equation-from-a-polynomial-surface-fit-to-a-positive-range?s_tid=srchtitle
instead of imposing constraints of f(xi)>=0, you imposes df/dx(xi) >= 0 for some xi "sufficiently" dense.
Deepa Maheshvare
Deepa Maheshvare 2022-9-6
@Bruno Luong, Thank you. I am not sure if I understand how to get k, l and C.
% cofficients of 2D polynomial 3d order
k = [0 0 1 0 1 2 0 1 2 3];
l = [0 1 0 2 1 0 3 2 1 0];
C = [xn.^k.*yn.^l]; % please no comment about my use of bracket here
% d = z;
Also could you please help with including
"instead of imposing constraints of f(xi)>=0, you imposes df/dx(xi) >= 0 for some xi "sufficiently" dense."
in the below code?
% ref: https://in.mathworks.com/matlabcentral/answers/1794985-how-to-constrain-the-resulting-equation-from-a-polynomial-surface-fit-to-a-positive-range#answer_1042065
x = [0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 20.85, 21.]
y = [2.83811035, 2.81541896, 3.14311655, 3.22373554, 3.43033456, 3.50433385, 3.66794514, 3.462296, 3.59480959,3.56250726, 3.6209845, 3.63034523, 3.68238915, 3.69096892, 3.75560395, 3.83545191, 3.90419498]
% Stuff needed to normalize the data for better inversion
[xmin, xmax] = bounds(x);
[ymin, ymax] = bounds(y);
xynfun = @(x,y)deal((x(:)-xmin)/(xmax-xmin),(y(:)-ymin)/(ymax-ymin));
[xn,yn]=xynfun(x,y);
% cofficients of 2D polynomial 3d order
k = [0 0 1 0 1 2 0 1 2 3];
l = [0 1 0 2 1 0 3 2 1 0];
C = [xn.^k.*yn.^l]; % please no comment about my use of bracket here
% d = z;
% Constraint positive of 3 x 3 points in the recatagular domain to be positive,
% it should be enough
[XNC,YNC] = meshgrid(linspace(0,1,3),linspace(0,1,3));
A = -[XNC(:).^k.*YNC(:).^l]; % please no comment ...
b = 0+zeros(size(A,1),1); % A*P<=0 means polynomial at (xnc,ync)>=0
P = lsqlin(C,d,A,b);
Bruno Luong
Bruno Luong 2022-9-6
编辑：Bruno Luong 2022-9-6
Ran in:
Ran in:
Monotonic polynomial
x = [0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 20.85, 21.];
y = [2.83811035, 2.81541896, 3.14311655, 3.22373554, 3.43033456, 3.50433385, 3.66794514, 3.462296, 3.59480959,3.56250726, 3.6209845, 3.63034523, 3.68238915, 3.69096892, 3.75560395, 3.83545191, 3.90419498];
% Stuff needed to normalize the data for better inversion
[xmin, xmax] = bounds(x);
xnfun = @(x)(x(:)-xmin)/(xmax-xmin);
xn=xnfun(x);
% cofficients of 2D polynomial 3d order
k = 0:7;
C = [xn.^k]; % please no comment about my use of bracket here
d = y;
% Constraint positive of 3 x 3 points in the recatagular domain to be positive,
% it should be enough
XNC = linspace(0,1,41);
A = -[k.*XNC(:).^(k-1)]; % please no comment ...
A(:,k==0)=0;
b = 0+zeros(size(A,1),1); % A*P<=0 means polynomial at (xnc,ync)>=0
P = lsqlin(C,d,A,b);
Minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance.
% Graphical check
% Create a grided model surface
xi=linspace(xmin,xmax,201);
Xin=xnfun(xi);
Yi=[Xin.^k]*P; % please no comment about my use of bracket here
close all
plot(xi,Yi);
hold on
plot(x,y,'or')
xlabel('x')
ylabel('y')

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Polynomials 的更多信息

产品

版本

R2021b

标签

Community Treasure Hunt

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

Start Hunting!

Translated by