Fitting a monotonically increasing spline function

40 views (last 30 days)
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.

Accepted Answer

Bruno Luong
Bruno Luong on 5 Sep 2022
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 Comments
Deepa Maheshvare
Deepa Maheshvare on 6 Sep 2022
@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 on 6 Sep 2022
Edited: Bruno Luong on 6 Sep 2022
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')

Sign in to comment.

More Answers (0)

Categories

Find more on Polynomials in Help Center and File Exchange

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!