lsqcurvefit Surface Fitting Troubleshooting

I am trying to write a script that fits any given surface using 2D polynomials. I am following the answer given in:

2D data fitting - Surface - MATLAB Answers - MATLAB Central

I am testing this by creating a surface using known polynomial coefficients, and then I can see if my script comes up with the correct coefficients through fitting it. Right now I'm using 9 coefficients. It seems my fitting works for the first three coefficients, but it ignores coefficients four through nine, and I cannot figure out why.

%% Setup
x = linspace(-100,100);
y = linspace(-100,100);
[X, Y] = meshgrid(x,y);
XY(:,:,1) = X; XY(:,:,2) = Y; 

C = [1 0 1 1 0 0 0 0 0]; %coefficients solution

Degree = 3;
Z = zeros(size(X)); %create original surface
idx = 1;
%define Z 
for j = 1:Degree
    for k = 0:j
        Z = Z + C(idx) * (X.^(j-k)) .* (Y.^k);
        idx = idx+1;
    end
end

%% Surface fitting 
% (from https://www.mathworks.com/matlabcentral/answers/119001-2d-data-fitting-surface)

poly_fun = @(c,XY) c(1)*XY(:,:,1).^1 + c(2)*XY(:,:,2).^1 + c(3)*XY(:,:,1).^2 + c(4)*XY(:,:,1).^1*XY(:,:,2).^1 + c(5)*XY(:,:,2).^2 + c(6)*XY(:,:,1).^3 + c(7)*XY(:,:,1).^2*XY(:,:,2).^1 + c(8)*XY(:,:,1).^1*XY(:,:,2).^2 + c(9)*XY(:,:,2).^3;
Cfit = lsqcurvefit(poly_fun, zeros(1,length(C)), XY, Z);

Zfit = zeros(size(X)); %Create fitted surface
idx = 1;
for j = 1:Degree
    for k = 0:j
        Zfit = Zfit + Cfit(idx) * (X.^(j-k)) .* (Y.^k);
        idx = idx+1;
    end
end

%% Plotting
figure
surf(Z,'LineStyle','none')
title('Original Surface')

figure
surf(Zfit,'LineStyle','none')
title('Polynomial Fit Surface')

Answer 1

[Torsten]

There is a problem with your "poly_fun". You have to use

poly_fun = @(c,XY) c(1)*XY(:,:,1).^1 + c(2)*XY(:,:,2).^1 + c(3)*XY(:,:,1).^2 + c(4)*XY(:,:,1).^1.*XY(:,:,2).^1 + c(5)*XY(:,:,2).^2 + c(6)*XY(:,:,1).^3 + c(7)*XY(:,:,1).^2.*XY(:,:,2).^1 + c(8)*XY(:,:,1).^1.*XY(:,:,2).^2 + c(9)*XY(:,:,2).^3;

instead of

poly_fun = @(c,XY) c(1)*XY(:,:,1).^1 + c(2)*XY(:,:,2).^1 + c(3)*XY(:,:,1).^2 + c(4)*XY(:,:,1).^1*XY(:,:,2).^1 + c(5)*XY(:,:,2).^2 + c(6)*XY(:,:,1).^3 + c(7)*XY(:,:,1).^2*XY(:,:,2).^1 + c(8)*XY(:,:,1).^1*XY(:,:,2).^2 + c(9)*XY(:,:,2).^3;

But this is a linear fitting problem in the coefficients c(1),...,c(9). Thus a nonlinear fitting tool like "lsqcurvefit" is not needed - you can simply use backslash (see the last 5 lines of the following code).

%% Setup
x = linspace(-100,100);
y = linspace(-100,100);
[X, Y] = meshgrid(x,y);
XY(:,:,1) = X; XY(:,:,2) = Y; 
C = [1 0 1 1 0 0 0 0 0]; %coefficients solution
Degree = 3;
Z = zeros(size(X)); %create original surface
idx = 1;
%define Z 
for j = 1:Degree
    for k = 0:j
        Z = Z + C(idx) * (X.^(j-k)) .* (Y.^k);
        idx = idx+1;
    end
end
%% Surface fitting 
% (from https://www.mathworks.com/matlabcentral/answers/119001-2d-data-fitting-surface)
poly_fun = @(c,XY) c(1)*XY(:,:,1).^1 + c(2)*XY(:,:,2).^1 + c(3)*XY(:,:,1).^2 + c(4)*XY(:,:,1).^1.*XY(:,:,2).^1 + c(5)*XY(:,:,2).^2 + c(6)*XY(:,:,1).^3 + c(7)*XY(:,:,1).^2.*XY(:,:,2).^1 + c(8)*XY(:,:,1).^1.*XY(:,:,2).^2 + c(9)*XY(:,:,2).^3;
Cfit = lsqcurvefit(poly_fun, zeros(1,length(C)), XY, Z);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

Zfit = zeros(size(X)); %Create fitted surface
idx = 1;
for j = 1:Degree
    for k = 0:j
        Zfit = Zfit + Cfit(idx) * (X.^(j-k)) .* (Y.^k);
        idx = idx+1;
    end
end
%% Plotting
figure
surf(Z,'LineStyle','none')
title('Original Surface')

●

figure
surf(Zfit,'LineStyle','none')
title('Polynomial Fit Surface')

X = X(:);
Y = Y(:);
Z = Z(:);
M = [X,Y,X.^2,X.*Y,Y.^2,X.^3,X.^2.*Y,X.*Y.^2,Y.^3];
c = M\Z

c = 9×1

1.0000 -0.0000 1.0000 1.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>