Constraining fit to two-variable function
Show older comments
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
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
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?
This seems to be inherently satisifed with the following constraints on the range occupied by your data, at least, and still seems to fit your data very well.
knownCoeffs= {'C00','C10','C20','C30','C40', 'C01','C02','C03','C04' };
knownVals=num2cell([1, zeros(1,8) ]);
If I want to make this a 5th order polynomial, I can just change ... and then set my coeffs that I need set correct?
In that case you would also need C50=0 in order to preserve the property that H=1 at eta=0.
How can I also extend it to eta/eta_c = 1?
With the same method. Just add points at eta/eta_c = 1 to the list of input points.
Benjamin
on 25 Mar 2019
Matt J
on 25 Mar 2019
You need to add corresponding points to extraR as well.
Matt J
on 26 Mar 2019
Yes.
Why does it become garbage when I copy the coeffs, the equation, and plot g(r) for a given eta value?
Because, you copied them only to the number of decimal places displayed on the screen, I assume. Also, because 6th order is a pretty high polynomial order, and so is getting to be somewhat unstable.
Edit: Ahh, maybe this: coeff=coeffvalues(fobj)
That is definitely what you would do instead of copying off the screen. However, a simpler way to restrict fobj to a specifc eta is to write an anonymous function
eta0=0.4/eta_c;
f_restricted = @(r) fobj(r,repelem(eta0,size(r)))
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?
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
More Answers (0)
Categories
Find more on Logical in Help Center and File Exchange
