Constraining fit to two-variable function

Question

Benjamin 2019-3-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/452213-constraining-fit-to-two-variable-function

编辑： Matt J 2019-3-28

采纳的回答： Matt J

example.xlsx

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2019-3-25

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/452213-constraining-fit-to-two-variable-function#answer_367285

编辑：Matt J 2019-3-25

在 MATLAB Online 中打开

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 个评论
显示 13更早的评论隐藏 13更早的评论

Benjamin 2019-3-26

编辑：Benjamin 2019-3-26

在 MATLAB Online 中打开

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

Benjamin 2019-3-27

编辑：Benjamin 2019-3-27

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 2019-3-28

编辑：Matt J 2019-3-28

在 MATLAB Online 中打开

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

请先登录，再进行评论。

Constraining fit to two-variable function

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

15 个评论
显示 13更早的评论隐藏 13更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Constraining fit to two-variable function

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

15 个评论 显示 13更早的评论隐藏 13更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

15 个评论
显示 13更早的评论隐藏 13更早的评论