MATLAB: Finding points of zero slope from 3D function

curve fittingCurve Fitting ToolboxImage Processing ToolboxsymbolicSymbolic Math Toolbox

  • I don't think this is what you want, but if you just want to plot the maximum you can use contourslice, i.e.
...
zeroPoints =...
% start added code
max_val = max(max(eval(DoG3DFit(X(1:100:end,1:100:end),Y(1:100:end,1:100:end)))));
[xmesh,ymesh,zmesh] = meshgrid(-1000:1000, -1000:1000,[0 max_val]);
V = matDoG3DFit(xmesh, ymesh);
hold on
cs = contourslice(xmesh,ymesh,zmesh,V,[],[],max_val,[max_val max_val]);
cs(1).EdgeColor = 'red';
cs(1).LineWidth = 3;
% end added code
subplot(2,2,3) ...
  • But from your question it looks as if you want to be able to generate an analytical solution.Doesn't the original (2D) method still work? I think at the 'caldera' both dz/dx = 0 and dz/dy = 0.So if you find either of these derivatives, you can find the maximum/minimum.Offcourse these derivatives are a bit more complicated..
Hope this helps.

Best Answer

To find the solution, you search the x and y for which the partial derivatives are null. This is tricky for an automated solver because the solution is complicated. A single solution in the center of the caldera, and an infinite number of solutions along the rim.
Below is an example of how you can calculate the symbolic solutions (R2017a).
syms a1 a2 b1 b2 b3 b4 c x x0 y y0
DoG3D(a1, a2, b1, b2, b3, b4, c, x, x0, y, y0) = ...
    (a1 * exp(-1/2 * (((x - x0)/b1)^2 + ((y - y0)/b2)^2))) - ...   % Gaus 1
    (a2 * exp(-1/2 * (((x - x0)/b3)^2 + ((y - y0)/b4)^2))) + ...   % Gaus 2
    c;
% The constants below will differ, depending on a fit, but here are some
% convenient ones
%                      a1   a2   b1  b2   b3   b4   c   x  x0   y  y0
DoG3DFit(x, y) = DoG3D(-50, -25, 50, 150, 100, 300, 10, x, 220, y, 180);
assume(x, 'real')
assume(y, 'real')
% Calculate the partial derivatives
dfdx = simplify(diff(DoG3DFit(x, y), x));
dfdy = simplify(diff(DoG3DFit(x, y), y));
% Solve with respect to y (3 solutions)
soly = simplify(solve(dfdy == 0, y));
% Substitute
dfdy = simplify(subs(dfdy, y, soly));
dfdx = simplify(subs(dfdx, y, soly));
% Only solution 1 has dependency in x
solx = simplify(solve(dfdx(1) == 0, x));
% Point solutions
pointSolutions = [solx repmat(soly(1), numel(solx), 1)]
% Solutions 2 and 3 are functions of x
functionSolutions = soly(2:3)
%%Plot the solutions
f2 = matlabFunction(soly(2), 'Vars', x);
f3 = matlabFunction(soly(3), 'Vars', x);
X = linspace(-500, 500, 3000);
Y2 = f2(X);
i2 = arrayfun(@(x) imag(x) == 0, Y2);
Y3 = f3(X);
i3 = arrayfun(@(x) imag(x) == 0, Y3);
figure;
fsurf(DoG3DFit, [-1000, 1000])
hold on
for k=1:size(pointSolutions, 1)
    plot3(pointSolutions(k, 1), pointSolutions(k, 2), DoG3DFit(pointSolutions(k, 1), pointSolutions(k, 2)), 'or', 'MarkerSize', 10, 'LineWidth', 5);
end
plot3(X(i2), Y2(i2), DoG3DFit(X(i2), Y2(i2)), 'm', 'LineWidth', 5);
plot3(X(i3), Y3(i3), DoG3DFit(X(i3), Y3(i3)), 'm', 'LineWidth', 5);
hold off
xlabel('x');
ylabel('y')
zlabel('z');
Related Question