# Constraining fit to two-variable function

10 views (last 30 days)
Benjamin on 25 Mar 2019
Edited: Matt J on 28 Mar 2019
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
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]);

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);

Matt J on 27 Mar 2019
1) Is there an easy way to get all the 2D plots? I'm currently doing it manually. Is there a way to loop through and produce all 2d plots with this function plotted against the data?
for i=1:N
f_restricted = @(r) fobj(r, repelem(eta(i)/eta_c,size(r)) ) ;
fplot(f_restricted)
end
2) How would I extend the 3d plot so that the surface goes from r=1 to 3, instead of 1 to 2? Is it just:
Did you try it?
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 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