Constraining fit to two-variable function

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

Answer 1

[Matt J]

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]

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

[Benjamin]

If I want to make this a 6th order polynomial, I can just change

 [I,J]=ndgrid(0:4);

to

[I,J]=ndgrid(0:6); 

and then set my coeffs that I need set correct? I could set the extra coeffs like:

knownCoeffs=       {'C00','C10','C20','C30','C40','C50','C60',   'C01','C11','C21','C31','C41','C51','C61'};
knownVals=num2cell([ [1   ,   0 ,   0 ,   0 ,   0 , 0, 0],   [ 8  ,   -6 ,   0 ,  0.5 ,   0 , 0 , 0]*eta_c ]);

Also, this does allow the fit to extend to eta/eta_c = 0. How can I also extend it to eta/eta_c = 1?

[Matt J]

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]

When I change it to:

extraR=[min(r);max(r)]; extraEta=[0;0;1;1]; extraH=fobj(extraR,extraEta/eta_c);

I get that matrix dimensions must agree. I'm not sure how to add this region.

[Matt J]

You need to add corresponding points to extraR as well.

[Benjamin]

Is it just:

extraR=[min(r);max(r);min(r);max(r)]; extraEta=[0;0;1;1]; extraH=fobj(extraR,extraEta/eta_c);

[Matt J]

Yes.

[Benjamin]

Are you sure the coefficients are output properly? If I take that equation and plot g(r) as a function of r, for a given eta (2D plot), g(r) behaves very oddly, and nothing like the fit. Is there a mistake in how the coeffs are output? The g(r) produced from the equation looks very different than the data when I take the output from your script. Also, note that I changed it to a 6-order fit. The fit looks good. Why does it become garbage when I copy the coeffs, the equation, and plot g(r) for a given eta value?

Could you try it? I have something like below. The coeffs that come with running it, produce a g that does not match the data at all. I imagine there must be a mistake in how the coeffs are being written out?

Note that I let some of the C-values vary that we had previously set. Basically I now allow C01, C11, and C31 to vary. All other Ci1's are set to 0.

Any idea why this produces a completely different function?

Note that I changed eta/eta_c in the code to just rpf, and set rpf earlier. Even keeping the code the same as you had it, I get the same results. Maybe if you try it, it will become obvious what is wrong?

Here is the code I used to get the fit:

close all;
clear;
clc;
%% Load Data
 num=xlsread('C:\example.xlsx',1,'A2:H18110');
eta_c= 0.6452;
r = num(:,4);
eta = num(:,3);
H = num(:,8);

rpf=eta/eta_c;
%% Set-up for fit
[I,J]=ndgrid(0:6); 
Terms= compose('C%u%u*r^%u*rpf^%u', I(:),J(:),I(:),J(:)) ;
  Terms=strjoin(Terms,' + ');
independent={'r','rpf'};
dependent='H';

knownCoeffs=       {'C00','C10','C20','C30','C40','C50','C60',   'C21','C41','C51','C61'};
knownVals=num2cell([ [1 ,   0 ,   0 ,   0 ,   0 ,   0 ,   0 ],   [ 0 ,  0 ,    0 ,   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)
  extraR=[min(r);max(r);min(r);max(r)]; extraEta=[0;0;0.6452;0.6452]; extraH=fobj(extraR,extraEta/eta_c);
  
  Rp=[extraR;r]; Etap=[extraEta; eta]/eta_c; Hp=[extraH;H] ;
  



hp=plot(fobj,[Rp,Etap],Hp);
% 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 'g(r,\eta)'
 zlim([0 1.1]);
 ylim([0 1.1]);

and here is the code that I took the fit into and tried to plot 2D plots. What is wrong with this? It is such a good fit, but then the 2d plots look really bad. Why would this be? Note that g(r) should be between 0 and 1. The fit when I put into MATLAB gives me a function that goes way above 1. This examples goes up to over 20

clear; clc; close all;
C01	=	3.66
C02	=	17.15
C03	=	-8.963
C04	=	-0.6612
C05	=	1059
C06	=	3610
C11	=	-4.024
C12	=	-25.22
C13	=	-0.1194
C14	=	-0.6537
C15	=	-3043
C16	=	-1.82E+04
C22	=	0.7829
C23	=	17.31
C24	=	88.37
C25	=	3014;
C26	=	3.69E+04;
C31	=	0.3634;
C32	=	9.708;
C33	=	-0.2074;
C34	=	-136.2;
C35	=	-1054;
C36	=	-3.89E+04;
C42	=	-0.2209;
C43	=	-11.08;
C44	=	-0.2726;
C45	=	0.07973;
C46	=	2.26E+04;
C52	=	-2.763;
C53	=	0.3303;
C54	=	81.49;
C55	=	0.09734;
C56	=	-6845;
C62	=	0.6396;
C63	=	2.247;
C64	=	-31.14;
C65	=	24.55;
C66	=	847.7;



       C00 =           1;
       C10 =           0;
       C20 =           0;
       C30 =           0;
       C40 =           0;
       C50 =           0;
       C60 =           0;
       C21 =           0;
       C41 =           0;
       C51 =           0;
       C61 =           0;  

eta=0.4;
       eta_c=0.6452;
       rpf=eta/eta_c;
       count = 0;
       for r=1:0.01:2
       count = count + 1;
        g(count) = C00*r^0*rpf^0 + C10*r^1*rpf^0 + C20*r^2*rpf^0 + C30*r^3*rpf^0 + ...
                    + C40*r^4*rpf^0 + C50*r^5*rpf^0 + C60*r^6*rpf^0 + C01*r^0*rpf^1 + ...
                    + C11*r^1*rpf^1 + C21*r^2*rpf^1 + C31*r^3*rpf^1 + C41*r^4*rpf^1 + ...
                    + C51*r^5*rpf^1 + C61*r^6*rpf^1 + C02*r^0*rpf^2 + C12*r^1*rpf^2 + ...
                    + C22*r^2*rpf^2 + C32*r^3*rpf^2 + C42*r^4*rpf^2 + C52*r^5*rpf^2 + ...
                    + C62*r^6*rpf^2 + C03*r^0*rpf^3 + C13*r^1*rpf^3 + C23*r^2*rpf^3 + ...
                    + C33*r^3*rpf^3 + C43*r^4*rpf^3 + C53*r^5*rpf^3 + C63*r^6*rpf^3 + ...
                    + C04*r^0*rpf^4 + C14*r^1*rpf^4 + C24*r^2*rpf^4 + C34*r^3*rpf^4 + ...
                    + C44*r^4*rpf^4 + C54*r^5*rpf^4 + C64*r^6*rpf^4 + C05*r^0*rpf^5 + ...
                    + C15*r^1*rpf^5 + C25*r^2*rpf^5 + C35*r^3*rpf^5 + C45*r^4*rpf^5 + ...
                    + C55*r^5*rpf^5 + C65*r^6*rpf^5 + C06*r^0*rpf^6 + C16*r^1*rpf^6 + ...
                    + C26*r^2*rpf^6 + C36*r^3*rpf^6 + C46*r^4*rpf^6 + C56*r^5*rpf^6 + ...
                    + C66*r^6*rpf^6;
       end
       r=1:0.01:2;
       plot(r,g)
       grid on;
       xlabel('x');
       ylabel('g(r)');
       
       
       

[Matt J]

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.

[Benjamin]

How do I copy more decimals? Do you get a good 2d plot? I need to actually use this fit so I need it to be accurate. It would seem that if the surface is such a good fit to the data between 0 and 1, it would be weird that when implementing the equation, I would get a function between 1 and 20. Perhaps I just need more decimals, but I'm not sure how to extract more.

Edit: I also have similar issues with 5th order. I think I just need more decimals. How do I get more decimals?

Edit: Ahh, maybe this:

 coeff=coeffvalues(fobj)

[Matt J]

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

[Benjamin]

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?

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:

extraR=[min(r);3;min(r);3]; extraEta=[0;0;0.6452;0.6452];

[Matt J]

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]

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]

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