Adding constraints to lsq fitting

I have a code that is of a function g(X,Y). The below code fits my data very well. However, I want to add some constraints at the boundaries where my data is not. For instance, how can I set the coefficient such that when Y=0, then the value of my function is 1?

Can anyone help with this?

clear
clc
close all

num=xlsread('C:data.xlsx',1,'A2:F18110');

eta = num(:,3);
r = num(:,4);
g = num(:,6);

%Do the surface fit
options=optimoptions(@lsqcurvefit,'SpecifyObjectiveGradient',false,'MaxFunctionEvaluations',50000,'MaxIterations',1000);

xdata={r,eta};
params0=linspace(0.1, 0.1, 50);
[params, resnorm]=lsqcurvefit(@modelfun,params0,xdata,g,[],[],options)

  for i=1:50
      C(i)=params(i);
  end
  
xmin = 1; xmax = 2; dx = 0.01;
etamin=0; etamax=0.55; deta=0.01;

[Xgrid,etagrid]=ndgrid(xmin:dx:xmax,etamin:deta:etamax);
surf(Xgrid,etagrid,modelfun(params,{Xgrid,etagrid}))
zlim([0 8]);
hold on;

scatter3(r(:),eta(:),g(:),'MarkerEdgeColor','none',...
         'MarkerFaceColor','b','MarkerFaceAlpha',.5); hold off
xlabel 'x', ylabel '\eta', zlabel 'g(x,\eta)'


function [out,Jacobian]=modelfun(params,xdata)
  X=xdata{1};
  Y=xdata{2};


  for i=1:50
      C(i)=params(i);
  end
     
    A1 = -0.4194;
    A2 = 0.5812;
    A3 = 0.6439;
    A4 = 0.4730;
    eta_c = 0.6452;
  
    PV = 1 + 3*Y./(eta_c-Y)+ A1*(Y./eta_c) + 2*A2*(Y./eta_c).^2 + 3*A3*(Y./eta_c).^3 + 4*A4*(Y./eta_c).^4;
    rdf_contact = (PV - 1)./(4*Y);

    poly_guess  = polyVal2D(C,X-1,Y/eta_c,4,4);
    
    out =  (poly_guess.*rdf_contact);
    
  if nargout>1
     
      %Jacobian=[X(:), Y(:), ones(size(X(:)))];
  end
  
end

Answer 1

[Adam Danz]

"...I want to add some constraints at the boundaries... "

In your call to lsqcurvefit(), you're not using the upper and lower bounds options (5th and 6th inputs are empty). I'd start out by imposing boundaries.

"For instance, how can I set the coefficient such that when Y=0, then the value of my function is 1"

I don't follow this part. The variable "Y" doesn't exist in your sample code. Also, what do you mean "the value of my function"? Do you mean the output of your nonlinear function or do you mean the coefficient estimates produced by the fit?

[Benjamin]

Y does exist, it's Y=xdata{2}; So essentially Y=eta. How do I impose that my coefficient estimates produce a surface function that equals 1 as eta (Y) goes to 0?

Basically, if I have polynomial terms like a_ij*x^i*y^j, then the constraint of y->0 then my function =1, means the remaining terms must add up to 1. This means that C_00 = 1 and C_i0 = 0. How can I do this in my code?

Also, in the bounds part, how do I specify which variable is being bounded?

[Matt J]

Also, in the bounds part, how do I specify which variable is being bounded?

In the syntax,

x = lsqcurvefit(fun,x0,xdata,ydata,LB,UB,options)

the inputs LB and UB are lower and upper bounds on your parameters. If you set LB(i)=UB(i)=somevalue, then x(i) will be constrained to equal somevalue. So if, for example, you want to force your first coefficient to be zero and all the rest unconstrained, you could do

LB=-inf(size(x0));
UB=+inf(size(x0));

LB(1)=0; UB(1)=0 % forces first unknown parameter to be zero

x = lsqcurvefit(fun,x0,xdata,ydata,LB,UB,options)

[Benjamin]

I'm not sure this accomplishes what I want though. I need the surface produced by the fit to equal 1 for all values of x, when eta or in this case Y =0. This only sets the first coefficient to 0, but then the the whole function needs to equal 1. Also, I believe the coefficient would need to be 1, not zero.

How can I set it up so that the entire solution is 1 when eta (Y) =0?

[Matt J]

"...for example..."

You already know what coefficients you want to constraint and to what values. Just apply the technique we've outlined to them.

[Benjamin]

Ok, I have:

params0=linspace(0.1, 0.1, 50);

LB=-inf(size(params0));
UB=+inf(size(params0));

LB(1)=1; UB(1)=1; % forces first unknown parameter to be zero

[params, resnorm]=lsqcurvefit(@modelfun,params0,xdata,g,LB,UB,options)

This sets the first coefficient = 1. Does the first coefficient correspond to C00 though for the MATLAB function you showed me? Or is it for a higher order term? How can I set all the Ci0 coefficients to 0?

[Adam Danz]

@Benjamin, sorry, I see the "Y" variable now.

I'm having trouble following the thread - mostly because I don't know which function you're refering to ("...the whole function needs to equal 1").

The upper/lower bound inputs to lsqcurvefit() will restrain the coeficient estimate outputs. Are you saying that you'd like a set of coeficients that, when used as inputs to your modelfun() function, will result in an output of 1?

[UPDATED] in response to your updated comment above.

LB(1)=1; UB(1)=1; forces the bound to be 1 (your inline comments says 0) on the first coefficient. This coefficient is corresponds to the first parameter of your nonlinear function.

[Benjamin]

Yeah sorry for the confusion. Basically, when Y=0, my produced surface needs to be = 0. This surface is produced from lsqcurvefit algorithm. One way I can accomplish this I believe is to set the coefficient C00 = 1. Then the coefficients Ci0 = 0, since they need to sum to 0. I believe I have C00 set to 1 above. How can I set the Ci0 (Y=0) coefficients = 0?

However, I am not sure if the first coefficient even corresponds to C00...

Edit: Yes, you have the right idea. My first coefficeint (assuming that corresponds to C00) needs to be 1. The coefficients Ci0 need to be 0. I am using the script polyVal2D. Does the first coefficient correspond to C00, or is it the last coefficient?

[Adam Danz]

Couldn't you just restrain that as well?

LB(2)=0; UB(2)=0; % forces 2nd parameter to be zero

[Benjamin]

Wouldn't there be several that need to be restrained? Basically I have Cij. anytime j=0, then Ci0 = 0. There would be several of them right? How do I even know which coefficent corresponds to what part of the polynomial? What order are they in? First instance, the solution I currently have contians 26 coefficients. What if I want to change this later to another higher order polynomial? Is there a way to make it so whenever eta = 0 in the coefficients, then that coefficient is 0, except when x is also 0. I know this is complicated...maybe there's not an easier way...

[Matt J]

What order are they in?

I believe we covered that in a previous thread. The help doc for polyFit2D explains the ordering of the coefficients.

[Benjamin]

Yes, that is what I believe I am looking at. If it is this, http://poquitopicante.blogspot.com/2013/03/horners-method-in-2d.html , then changing the 1st coefficient to 1, is not correct right? I need to be changing the last coefficient.

[Adam Danz]

Assuming I've found the correct polyVal2D() function you're using,

% Polynomial coefficients are in the following order.
%
%   F(X,Y) = P_1 * X^N * Y^M + P_2 * X^{N-1} * Y^M + ... + P_{N+1} * Y^M + ...
%            P_{N+2} * X^N * Y^{M-1} + P_{N+3} * X^{N-1} * Y^{M-1} + ... + P_{2*(N+1)} * Y^{M-1} + ...
%            ...
%            P_{M*(N+1)+1} * X^N + P_{M*(N+1)+2} * X^{N-1} + ... + P_{(N+1)*(M+1)}

[Benjamin]

So then it is the last coefficient that needs to be set to 1 right?

Answer 2

[Matt J]

Here is an adaptation of the algebraic solution that I presented in your previous thread. I don't bother with the lsqcurvefit alternative, since as I showed you, it is both slower and less accurate. Below are also surface plots of the fitted surface, both in the neighborhood of the data points and near eta=0 where you can see that everything converges to the target value parameter, g0, as eta-->0. In this example, I have set g0=1.3.

clear
clc
close all


%% Data Set-up
num=xlsread('example.xlsx',1,'A2:F18110');
Np=4; %polynomial order
g0=1.3; %target value at eta=0

g = num(:,6);
r = num(:,4);
eta = num(:,3);

otherStuff.r = r;
otherStuff.eta = eta;
otherStuff.g = g;
otherStuff.eta_c = 0.6452;
otherStuff.Np=Np;
otherStuff.A1 = -0.4194;
otherStuff.A2 = 0.5812;
otherStuff.A3 = 0.6439;
otherStuff.A4 = 0.4730;

%% Fit using matrix algebra

A0=func2mat(@(p) modelfun(p,otherStuff), ones(Np+1,Np+1));

 q=4*g0./((3+otherStuff.A1)/otherStuff.eta_c); %weight on limiting target value as eta-->0
 
 b=g(:)-q*A0(:,1);
  A=A0;
 A(:,1:Np+1)=[];

 coeffs =[zeros(Np+1,1); A\b]; %fitted coefficients
 coeffs(1)=q;

figure(1); showFit(coeffs,otherStuff);

%% Also check surface near eta=0

  [Rg,Eg]=ndgrid(unique(r),linspace(1e-8,0.01,100));
  
 differentStuff=otherStuff;
 differentStuff.r = Rg(:);
 differentStuff.eta = Eg(:);
 differentStuff.g=nan([numel(Rg),1]);

figure(2); showFit(coeffs,differentStuff);


function ghat=modelfun(C,otherStuff)


        r = otherStuff.r;             
        eta = otherStuff.eta;                  
        eta_c = otherStuff.eta_c;     
        Np = otherStuff.Np;           
        A1 = otherStuff.A1;           
        A2 = otherStuff.A2;           
        A3 = otherStuff.A3;           
        A4 = otherStuff.A4;   

  
    PV = 1 + 3*eta./(eta_c-eta)+ A1*(eta./eta_c) + 2*A2*(eta./eta_c).^2 +...
         3*A3*(eta./eta_c).^3 + 4*A4*(eta./eta_c).^4;
     
    rdf_contact = (PV - 1)./(4*eta);

    Cflip=rot90(reshape(C,Np+1,Np+1),2);
    
    poly_guess  = polyVal2D(Cflip,r-1,eta/eta_c,Np,Np);
    
    ghat = (poly_guess.*rdf_contact);
    

end

function showFit(coeffs,otherStuff)
  
     r = otherStuff.r(:);        
    eta = otherStuff.eta(:);         
    g = otherStuff.g(:);           
    ghat=modelfun(coeffs,otherStuff);
    
    tri = delaunay(r,eta);
    %% Plot it with TRISURF
    h=trisurf(tri, r, eta, ghat);
    h.EdgeColor='b';
    h.FaceColor='b';
    axis vis3d
    
    hold on;
    scatter3(r,eta,g,'MarkerEdgeColor','none',...
             'MarkerFaceColor','r','MarkerFaceAlpha',.05); 
    xlabel 'r', ylabel '\eta', zlabel 'g(r,\eta)'
    hold off
end

[Benjamin Cowen]

Wow this is impressive! I’m not on my computer now so it’s hard to test it out myself, but it looks good! I’ll have to go through line by line to see what you did because it looks pretty different than what I had before. Could you explain why this method is better than lsq method? How is it even solving for the coefficients here if it is not iterative?

[Matt J]

Your model function is linear (in the coefficients). Therefore, there exists a matrix A and vector b such that the coefficients c is the solution to,

This has an analytic solution that is computed in Matlab via c=A\b. Note that this only works because your model function is linear. If you ever tweak the model making it non-linear you would have to go back to lsqcurvefit.

[Benjamin]

1) What exactly is the g0=1.3? I need my function to equal 1 at eta = 0. Note this does work since the 4eta divided actually cancels with part of the numerator. It is not clear to me exactly how the weighting is done. Can't we just set the coefficient to 1 and the other coefficients to 0 that are of the form I mentioned? Why is it necessary to have MATLAB solve for it? I am not sure where this equation comes from:

q=4*g0./((3+otherStuff.A1)/otherStuff.eta_c);

Can't we just set the C00 coefficeint to 1 and the Ci0 coefficients to 0? I think just setting the coefficients would be preferable, if possible. I would prefer to remove any type of weighting function and just set the coeffs. I also have other coefficients that I need to set. I see that you have q = all that stuff, but why can't we just set q =1 ? Then coeffs(1) = q. This seems like it does what I want right? Then can I just set the other coeffs? I'm still not sure what order they are in. Aren't the coeffs in reverse order?

If you could tell me explictly, what term in the C matrix is C00, C11, C12, etc, that would super helpfu. As it stands, I have no idea which is which.

2) How would I add in more constraints? The main last constraint that I need to add is that any terms that are eta^1 * x^i need to be fixed. If i = x and j = eta, then C10 = 8, C11 = -6, C13 = 1/2 and all other C1i = 0. Can this be done?

So to be clear. I need (where x = i, eta = j)

Cij format

C00 = 1;

Ci0's (other than the one listed above) = 0

C01 = 8

C11 = -6

C31 = 1/2

Ci1's (that are not listed above) = 0

3) This code creates an array called coeffs. I can see that the first entry is 1. This implies that the first entry is the C00 coefficient. According to http://poquitopicante.blogspot.com/2013/03/horners-method-in-2d.html isn't the last term supposed to be this coefficient?

4) In the plotted figure 1, how can I make the surface extend to eta = 0? Currently, it stops where the data does.

5) In

A0=func2mat(@(p) modelfun(p,otherStuff), ones(Np+1,Np+1))

what is p? It looks like its the coeffs, but where is this set?

I am especially confused over the coeffs.==

I just want to reiterate that I am very greatful for all of your help!

If the coeffs go like C00, C10, C20, C30, C40, C01, C11, then I added this code:

 coeffs(1)=1;
 coeffs(2)=0;
 coeffs(3)=0;
 coeffs(4)=0;
 coeffs(5)=0;
 coeffs(6)=8;
 coeffs(7)=-6;
 coeffs(9)=0.5;

But I am not sure if this is right. Now the surface fit seems to be way off.

[Matt J]

If you could tell me explictly, what term in the C matrix is C00, C11, C12, etc, that would super helpfu. As it stands, I have no idea which is which.

In modelfun, I have flipped the input to polyVal2D to make the ordering of the coefficients is more intuitive. In my implementation C(i,j) corresponds to the mononomial term r.^(i-1).*eta.^(j-1), for i,j=1,...,Np+1. So, it is esentially set up now as you've suggested. By making C(1,1)=constant and C(2:Np+1,1)=0, the surface is made to converge to a common value as eta-->0 for all values of r.

1) What exactly is the g0=1.3? I need my function to equal 1 at eta = 0.

g0 is an adjustable parameter. If you want the function to approach 1 as eta-->0, then set g0=1. Note that your surface cannot actually be evaluated at eta=0 (as currently coded) because the weighting term rdf_contact is infinite at eta=0. I think that you can fix that by re-implenting rdf_contact as below. This gets rid of all divisions by eta.

    rdf_contact = 3./(eta_c-eta)+ A1*(1./eta_c) + 2*A2*eta./eta_c.^2 +...
         3*A3*eta.^2./eta_c.^3 + 4*A4*eta.^3/eta_c.^4;
     
    Cflip=rot90(reshape(C,Np+1,Np+1),2);
    
    poly_guess  = polyVal2D(Cflip,r-1,eta/eta_c,Np,Np);
    
    ghat = (poly_guess.*rdf_contact);

Can't we just set the C00 coefficeint to 1 and the Ci0 coefficients to 0?

The first coefficient cannot be set to 1 because it must account for additional weighting contributed by rdf_contact.

2) How do I add more constraints...So to be clear. I need (where x = i, eta = j)

Didn't follow that. The coefficients themselves are not functions of r and eta. They are coefficients of things that are functions of r and theta.

3) isn't the last term supposed to be this coefficient?

As I mentioned, I flipped the coefficients. That's what this line does

        Cflip=rot90(reshape(C,Np+1,Np+1),2);

4) In the plotted figure 1, how can I make the surface extend to eta = 0? Currently, it stops where the data does.

This should do it. Now you can pass additional points [R,Eta] to showFit, which you would like added to the surface:

function showFit(coeffs,otherStuff, R,Eta)
  
     r = otherStuff.r(:);        
    eta = otherStuff.eta(:);         
    g = otherStuff.g(:); 
    
    otherStuff.r=[r;R(:)]; 
    otherStuff.eta=[eta;Eta(:)]; 
    ghat=modelfun(coeffs,otherStuff);
     
    
    tri = delaunay(otherStuff.r,otherStuff.eta);
    %% Plot it with TRISURF
    h=trisurf(tri, otherStuff.r, otherStuff.eta, ghat);
    h.EdgeColor='b';
    h.FaceColor='b';
    axis vis3d
    
    hold on;
    scatter3(r,eta,g,'MarkerEdgeColor','none',...
             'MarkerFaceColor','r','MarkerFaceAlpha',.05); 
    xlabel 'r', ylabel '\eta', zlabel 'g(r,\eta)'
    hold off
end

[Benjamin]

Can you remove the weighting function, and then I can just set my coeffs to what they need to be? I think the coeffs should look like this below. The weird thing is now the surface fit is really bad, so I must be missing something? In modefun(p ...), what is p? It looks like it should be the coeffs, but where is this stated?

 coeffs(1)=1;
 coeffs(2)=0;
 coeffs(3)=0;
 coeffs(4)=0;
 coeffs(5)=0;
 coeffs(6)=8;
 coeffs(7)=-6;
 coeffs(9)=0.5;

How do i pass more values or r and eta? I don't want to pass more data, I just want the surface to extend further out.

Do I add like r = [1 2] and eta=[0 1]? This gives me error of too many input arguments.

One idea too: I take my data and divide out rdf_contact, then I just have my polynomial in the fit. Can we try that?

I have no idea where this came from:

 q=4.*g0./((3+otherStuff.A1)/otherStuff.eta_c);

I took my data and divided out the rdf_contact. Now the function fit can completely remove the rdf_contact. I have attached new excel file. Just change g is equal to data in column H.

Edit: updated excel file

[Matt J]

Can you remove the weighting function, and then I can just set my coeffs to what they need to be?

It's your model. I assume rdf_contact has a reason to be there. As you'll recall, that was the whole reason you couldn't simply do this with straight up polynomial fitting with the Curve Fitting Toolbox.

Do I add like r = [1 2] and eta=[0 1]?

You add whatever additional points you want included in the surface plot.

I have no idea where this came from:

By setting C(1,1) to some constant q, as you've been proposing, and C(2:Np+1,1)=0 and letting eta-->0 you will find that the surface approaches the limiting value

LimitingValue=(3+otherStuff.A1)/otherStuff.eta_c/4 * q;

Therefore, setting

q=4.*g0./((3+otherStuff.A1)/otherStuff.eta_c));

leads to

LimitingValue = g0

[Benjamin]

So, I've changed the excel file. I now have already divided out rdf_contact. Here is the new excel file. The new g is column H. It already had divided out rdf_contact so I am just fitting the polynomial.

Any chance you could modify the code without the weighting function and with the new data? Note that this means rdf_contact and all the coeffs (A1, A2, ..) with it can be removed.

I apologize if I should have done this from the beginning.

When I use cftool on (r,eta,g), I feel like the fits don't look as good...It actually looks a lot worse than what you were able to do in your code. Unless I am missing something. I also don't feel like I have as much control over changing things.

Could I also add data points at eta =0 , that g=1 for all x? That might help the fit too...I keep trying wtih cftool, and the fits dont look as good as in your matrix algebra code. Any chance you could change the code with the new things? I still can use the old code for when I have rdf_contact. But I want to try it without.

So to sum up, new excel file, column H has new g data with rdf_contact already accounted for. Rdf_contact can be removed from matlab code. And I want to set:

So to be clear. I need (where x = i, eta = j)

Cij format

C00 = 1;

Ci0's (other than the one listed above) = 0

C01 = 8/eta_c

C11 = -6/eta_c

C31 = 1/2/eta_c

Ci1's (that are not listed above) = 0

So it would be good to remove weighting function I think and just set these coeffs. Then I should be good to go!

Edit: I keep trying to get good fits in Cftool... I think your matrix algebra is better. Any chance you could update the code to ignore the rdf_contact, just read in column h for g(x,eta), and set the above coefficients?

Edit: When I change your code, my surface fit is completely messed up. I'm not sure if I'm not putting the coefficients in correctly?

Does coeff(9) correspond to C31?

Also with these coeffs, r-1 needs to be just r.

Sorry for all the comments...

Also, you pass p to model fun. What is p? It's not defined anywhere that I can tell?

Answer 3

[Matt J]

Any chance you could modify the code without the weighting function and with the new data? Note that this means rdf_contact and all the coeffs (A1, A2, ..) with it can be removed.

It's virtually trivial now, with the Curve Fitting Toolbox:

%% Load Data
num=xlsread('example2.xlsx',1,'A2:F18110');
eta_c= 0.6452;

r = num(:,4);
eta = num(:,3);
H = num(:,5);

%% Set-up for fit

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

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],H,ft,'problem',knownVals) ;

hp=plot(fobj,[r,eta],H);
 set(hp(1),'FaceAlpha',0.5);
 set(hp(2),'MarkerFaceColor','r');
 xlabel 'r',
 ylabel '\eta'
 zlabel 'H(r,\eta)'

[Benjamin]

H should be column 8, not 5 I believe. column 8 is g without rdf_contact. Can the surface be exteneded to eta=0? How do I do that? I can change the limits, but how do I extend the surface to eta = 0?

Also, if I plot this, you can see the the curves close to eta =0.1 don't really match the fit that well. Any way to improve that?

Also, it should be a fit to eta/eta_c, not eta. eta_c=0.6452

When I try to change this:

fobj = fit([r,eta],H,ft,'problem',knownVals)

to


fobj = fit([r,eta/eta_c],H,ft,'problem',knownVals)

the resulting fit looks terrible.

In the old code, we had:

poly_guess  = polyVal2D(Cflip,r-1,eta/eta_c,Np,Np);

How do I add the /eta_c here? When I try it, the surface gets all messed up. Does it for you?

Also, those few coeffs that are divided by eta_c, should actually be multiplied by eta_c.

Thanks for this!

[Benjamin Cowen]

I still cannot fully get this figured out, but I feel like I am close. During the fit, is there a way to fit the function to x and eta/eta_c. Currently it is fit to x and eta. I changed in my code a few things. For instance, the coefficients that you have divided by eta_c, should actually be multiplied. Also, you have H equal to the 5th column, but it should be the 8th. When I make all these changes, the surface fit is pretty far off, but maybe I'm doing something incorrectly. Any chance you could take a look? I edited the code to fix the coefficients and change H to the 8th column. The updated code is below. How can I make the fit be to x and eta/eta_c? Mabye the curve looked off because the eta axis needs to be changed to eta/eta_c? Lastly, it would be nice if the code could be updated so that I can see the surface as eta->0.

Edit: I can't actually test it right now since I'm not on computer. I think maybe I realized that the axes was wrong. Since I fit to eta/eta_c, then the y-axis should be eta/eta_c I would think. Can you test this? Does the curve fit well with the code below?

I changed the things I mentioned and these 2 lines (changed eta to eta/eta_c).

fobj = fit([r,eta/eta_c],H,ft,'problem',knownVals) ;
hp=plot(fobj,[r,eta/eta_c],H);

here is the whole code.

%% Load Data
num=xlsread('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:4); 
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'
 zlabel 'H(r,\eta)'

[Matt J]

The fit that I get looks pretty good to me...

[Benjamin]

If you notice, the surface still does not extend to eta=0. How do I get it to do that? Currently, the surface is only plotted over the data. I would like the surface to go all the way to eta/eta_c = 0.

Edit: I also created a new question that is related to this in another thread.