Constraining fit to two-variable function

1 view (last 30 days)
Hello,
I have the following code, which produces a surface of my data. (Example xlsx file attached). Certain coeffs are already constrained. My question is, how do I constrain the fit so that H(r,eta) does not go above 1? Also, even though I do not have data past eta/eta_c of around 0.15, how do I extend my fitted surface to eta/eta_c = 0? At 0, the function should equal 1 since the C00 coefficient is set to 1.
close all
%% Load Data
num=xlsread('C:example.xlsx',1,'A2:H18110');
eta_c= 0.6452;
r = num(:,4);
eta = num(:,3);
H = num(:,8);
%% Set-up for fit
[I,J]=ndgrid(0:5);
Terms= compose('C%u%u*r^%u*eta^%u', I(:),J(:),I(:),J(:)) ;
Terms=strjoin(Terms,' + ');
independent={'r','eta'};
dependent='H';
knownCoeffs= {'C00','C10','C20','C30','C40', 'C01','C11','C21','C31','C41'};
knownVals=num2cell([ [1 , 0 , 0 , 0 , 0 ], [ 8 , -6 , 0 , 0.5 , 0 ]*eta_c ]);
allCoeffs=compose('C%u%u',I(:),J(:));
[unknownCoeffs,include]=setdiff(allCoeffs,knownCoeffs);
ft=fittype(Terms,'independent',independent, 'dependent', dependent,...
'coefficients', unknownCoeffs,'problem', knownCoeffs);
%% Fit and display
fobj = fit([r,eta/eta_c],H,ft,'problem',knownVals);
hp=plot(fobj,[r,eta/eta_c],H);
set(hp(1),'FaceAlpha',0.5);
set(hp(2),'MarkerFaceColor','r');
xlabel 'r',
ylabel '\eta/\eta_c'
zlabel 'H(r,\eta)'
zlim([0 1.1]);
ylim([0 0.55]);

Accepted Answer

Matt J
Matt J on 25 Mar 2019
Edited: Matt J on 25 Mar 2019
My question is, how do I constrain the fit so that H(r,eta) does not go above 1?
Polynomials cannot be bounded everywhere. You would at the very least have to decide on a particular region where this bound would apply.
At 0, the function should equal 1 since the C00 coefficient is set to 1.
For that, you need to fix this line:
[I,J]=ndgrid(0:4);
how do I extend my fitted surface to eta/eta_c = 0
Just pass extra points to the plot command,
extraR=[min(r);max(r)]; extraEta=[0;0]; extraH=fobj(extraR,extraEta/eta_c);
Rp=[extraR;r]; Etap=[extraEta; eta]/eta_c; Hp=[extraH;H] ;
hp=plot(fobj,[Rp,Etap],Hp);
  15 Comments
Benjamin
Benjamin on 27 Mar 2019
Edited: Benjamin on 27 Mar 2019
How do I add this to the code? When I replace the old plot with this one, I get Undefined function or variable N. If I replace N with a number, I get a blank plot. Is there a way to add this to the old code, where I get the 3d plot and all the 2d plots?
Matt J
Matt J on 28 Mar 2019
Edited: Matt J on 28 Mar 2019
Hmmm. Simpler than I thought.
Etas=0.2:0.2:0.8 ;
for i=1:numel(Etas)
f_restricted = @(r) fobj(r(:).', Etas(i)/eta_c ) ;
figure(i), fplot(f_restricted)
end

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!