Newtons method with two variables

So I'm trying to implement Newtons method to find the coordinates of an extremevalue for a function h, whose first and second partial derivates I've named f and g respectively. My code does not return the x and y-coordinates, and I believe the fault lies in my if-statement. Here's the code:

if true
  [x,y] = meshgrid (-3:0.1:3,-3:0.1:3);
h = @(x,y) ((x.^3.*y) + ((5*x.^2).*(y.^2)))./(exp(1).^((x.^2) + 3*y.^4));
f = @(x,y) ((3*x.^2.*y)+(10*y.^2.*x)-((10*y.^2).*(x.^3))-(2*x.^4*y))./(exp(1).^((x.^2)+3*y.^4));
%h är funktionen given i uppgiften. f är den partiella derivatan av h med
%avseende på x. Denna omskriving görs på grund av förenklad notation.
g = @(x,y) (((x.^3)+(10*x.^2.*y))-((60*x.^2).*(y.^5))+((12*x.^3).*(y.^4)))./(exp(1).^((x.^2)+3*y.^4));
%h är funktionen given i uppgiften. g är den partiella derivatan av h med
%avseende på y. Denna omskriving görs på grund av förenklad notation.
%contour(x,y,f(x,y)) (contour och mesh skapar graferna i uppgift 2a. Samma
                    %kod används här som den i uppgift 2a som skapade 
                    %figurerna i labbrapporten.) 
%mesh(x,y,f(x,y))
f1 = @(x,y) ((4*x.^5.*y) + ((20*x.^4).*(y.^2))-(14*x.^3.*y)-((50*x.^2).*(y.^2))+(6*x.*y)+(10*y.^2))./(exp(1).^((x.^2)+3*y.^4));
f2 = @(x,y) ((((120*x.^3).*(y.^5))-(120*y.^5.*x))+(((24*x.^3).*(y.^4))-(36*y.^4.*x))+(20*y-(20*x.^2.*y))-2*x.^3+3*x)./(exp(1).^((x.^2)+3*y.^4));
g1 = @(x,y) ((((24*y.^4).*(x.^4))-(2*x.^4))+(((120*y.^5).*(x.^3))-((20*y.^2).*(x.^2))+(3*x.^2)-((36*y.^4).*(x.^2))-(120*y.^5.*x)+(10*y.^2)))./(exp(1).^((x.^2)+3*y.^4));
g2 = @(x,y) (((720*y.^8).*(x.^2))+((144*y.^7).*(x.^3))-(120*y.^5.*x)-((300*y.^4).*(x.^2))-((60*x.^3).*(y.*3))+(20*x.*y))./(exp(1).^((x.^2)+3*y.^4));
x(1) = 1.1;
y(1) = 0.6;
h = 0.0000001;
X = [f(x(i),y(i)) f2(x(i),y(i)); g(x(i),y(i)) g2(x(i),y(i))];
Y = [f1(x(i),y(i)) f(x(i),y(i)); g1(x(i),y(i)) g(x(i),y(i))];
J = [f1(x(i),y(i)) f2(x(i),y(i)); g1(x(i),y(i)) g2(x(i),y(i))];
for i=1
  x(i+1) = x(i) - (det(X)/det(J)); %([f(x(i),y(i)) f2(x(i),y(i)); g(x(i),y(i)) g2(x(i),y(i))] / [f1(x(i),y(i)) f2(x(i),y(i));g1(x(i),y(i)) g2(x(i),y(i))]);
  x(i) = x(i+1);
  y(i+1) = y(i) - (det(Y)/det(J));
  y(i) = y(i+1);
  if (f(x(i),y(i)) && g(x(i),y(i)) < h )
      valX = x(i);
      valY = y(i);
      return;
  end
end
figure()
hold on;
contour(x,y,f(x,y),[0 0], 'b')
contour(x,y,g(x,y),[0 0], 'r')
end

Note: the final four rows simply draw a figure that shows approximate critical points. Comments in swedish. Thanks!

Answer 1

[John D'Errico]

Why in the name of god and little green apples do you need to write your own code?

NEVER write your own code, when code written by professionals is available. This is especially true when you will write poor code.

ezsurf(h,[-3,3])

Rotate the plot around, and it will be obvious the maximum lies near (1,1). It will also be clear there will be multiple solutions. As well, it looks like there might be two solutions with the same maximum, mirror images across the origin.

fminsearch(@(xy) -h(xy(1),xy(2)),[1 1])
ans =
        1.062      0.61742

Or, you could have used the optimization toolbox, which will probably give you a few more correct digits for the time required.

format long g
fminunc(@(xy) -h(xy(1),xy(2)),[1 1])

Local minimum found.

Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.

<stopping criteria details>

ans =
        1.06206524800884         0.617438162285868

Next, is it true that the two solutions are the mirrored across the origin?

[xy1,fval] = fminsearch(@(xy) -h(xy(1),xy(2)),[-1 -1])
xy1 =
         -1.06203291314256         -0.61741902622501
fval =
        -0.604829999447586

[xy1,fval] = fminsearch(@(xy) -h(xy(1),xy(2)),[1 1])
xy1 =
          1.06203291314256          0.61741902622501
fval =
        -0.604829999447586

Looks like it.

Or, you could use the symbolic TB.

syms x y
H = ((x.^3.*y) + ((5*x.^2).*(y.^2)))./(exp(1).^((x.^2) + 3*y.^4));

Hgrad = gradient(H);

xysol = solve(Hgrad,x,y)
xysol = 
  struct with fields:

      x: [13×1 sym]
      y: [13×1 sym]

So 13 solutions, many of which are complex. Which one givex the max?

vpa([xysol.x,xysol.y,h(xysol.x,xysol.y)],8)
ans =
[                 0.90244248,                 -0.66672465,                 0.32319386]
[                 -1.0620653,                 -0.61743817,                    0.60483]
[                  1.4133748,                 -0.14100149,               -0.027034788]
[                 -1.4133748,                  0.14100149,               -0.027034788]
[                  1.0620653,                  0.61743817,                    0.60483]
[                -0.90244248,                  0.66672465,                 0.32319386]
[   1.0188998 + 0.076070383i, - 0.025239699 - 0.63521774i, - 0.44731804 - 0.14437493i]
[   1.0188998 - 0.076070383i, - 0.025239699 + 0.63521774i, - 0.44731804 + 0.14437493i]
[                  10.20644i,                 -2.0509926i,               1.5918468e23]
[                 -10.20644i,                  2.0509926i,               1.5918468e23]
[ - 1.0188998 + 0.076070383i,   0.025239699 - 0.63521774i, - 0.44731804 + 0.14437493i]
[ - 1.0188998 - 0.076070383i,   0.025239699 + 0.63521774i, - 0.44731804 - 0.14437493i]
[                          0,                           0,                          0]

Of the real solutions found, it looks like the two we found before are it.