Hessian matrix 3D plot help

I am writing a script to find/plot the maxs and mins of a 3D function using the steepest ascent method, golden search method and confirmed by the Hessian matrix, I have almost finished this question for a school project but am having trouble having with the script storing the appropriate max/min/saddle values to plot on the 3D surface. The script plots the 3D surface its self, but no max or min points. I have attached what I have created so far, I just need some pointers in allocating the max/min/saddle points so they will be plotted on the 3D surface.

clear;clc

Fun = @(q,r)1.8.^(-1.5*sqrt(q.^2+r.^2)).*cos(0.5*r).*sin(q);

Guess = [2,2; 2,-2; -2,2; -2,-2];                      %Give intial guesses

XRange = [-2 2];         YRange = [-2 2];

%------------Establish a matrix of four starting guesses based on the range-------------%

x = XRange(1,1) : 0.25 : XRange(1,2);                  %(-2, -1.75, -1.5,...,1.5, 1.75, 2)%
y = YRange(1,1) : 0.25 : YRange(1,2);

[X,Y] = meshgrid(x,y);                                 %Form a mesh of base points that the surf plot will use
Z = 1.8.^(-1.5*sqrt(X.^2+Y.^2)).*cos(0.5*Y).*sin(X);
deltax = 0.01;                                        % Establish delta's for central difference method
deltay = 0.01;

iMin = 1; iMax = 1; iSad = 1;                          %Keep the number of mins, maxs, and saddle points found
Min = zeros(1,5); Max = zeros(1,5); Saddle = zeros(1,5);
d = zeros(1,5); h_opt = zeros(1,5); xopt = zeros(1,5);
f1 = zeros(1,5); f2 = zeros(1,5);
x1 = zeros(1,5); x2 = zeros(1,5);
xU = zeros(1,5); xL = zeros(1,5);
GS_Error = zeros(1,5); SA_Error = zeros(1,5);
R = (5^.5-1)/2; 

%---------Steepest Ascent Method-------%
for w = 1:4                                            %Loop over the four intial guesses
    x_new = Guess(w,1); y_new = Guess(w,2);            %Location of intial guesses
    OptZ_new = Fun(x_new, y_new);                      %Funcion value at the point

      while SA_Error>.01
          x_old = x_new;
          y_old = y_new;
          OptZ_old = OptZ_new;
          InitialGuess=Fun(x_new,y_new);
          dfdx = (Fun(x_old+deltax,y_old)-Fun(x_old-deltax,y_old))/(2*deltax);
          dfdy = (Fun(x_old,y_old+deltay)-Fun(x_old,y_old-deltay))/(2*deltay);
          Funh = @(h)Fun(x_new+dfdx*h,y_new+dfdy*h);  

          %Begin Golden Section method

          while abs(GS_Error)>.01

                  d(i) = ((sqrt(5)-1)/2)*(xU(i)-xL(i));
                  x1(i)=xL(i)+d(i);
                  x2(i)=xU(i)-d(i);
                  f1(i)=Funh(x1(i));
                  f2(i)=Funh(x2(i));

                  if abs(f1(i))>abs(f2(i))
                      OptZ_high = x1;
                      x1 = x2;
                      x2 = OptZ_high - d(i);
                  elseif abs(f2(i))>abs(f1(i))
                      OptZ_low= x2;
                      x2 = x1;
                      x1 = OptZ_low +d(i);
                      GS_Error(i) = (1-R)*abs(xU(i)-xL(i))/xopt(i);
                      x_new=x_new+dfdx*h_opt;                         
                      y_new=y_new+dfdy*h_opt;
                      SA_Error = sqrt((x_new-x_old)^2+(y_new-y_old)^2);
                  end
                  i = i+1;
              end
          end
  end

%Use Hessian matrix to determine local minimum or maximum
d2fdx2 = (Fun(x_new+deltax,y_new)-2*Fun(x_new,y_new)+Fun(x_new-deltax,y_new))/(2*deltax)^2;
d2fdy2 = (Fun(x_new,y_new+deltay)-2*Fun(y_new,x_new)+Fun(x_new,y_new-deltay))/(2*deltay)^2;
d2fdxdy = (Fun(x_new+deltax,y_new+deltay)-Fun(x_new+deltax,y_new-deltay)-Fun(x_new-deltax,y_new+deltay)+Fun(x_new-deltax,y_new-deltay))/(4*deltax*deltay);
Hdet=d2fdx2*d2fdy2-d2fdxdy^2;

if abs(Hdet) > 0 && d2fdx2 > 0;
       Max = OptZ_new;
elseif abs(Hdet) > 0 && d2fdx2 < 0;
       Min = OptZ_new;
elseif abs(Hdet) <0;
    saddle = OptZ_new;
end

%Plot for found mins and maxs
hold on
surf(X,Y,Z)
plot3(Max(1,1),Max(1,2),Min(1,1),Min(1,2),Saddle(1,1),Saddle(1,2),'m*','MarkerSize',10)

Thank you in advance for any help or advice, I feel like it is something very simple I am just not seeing.

Answer 1

[James Tursa]

If you just need help with plotting the point, maybe something like this:

if abs(Hdet) > 0 && d2fdx2 > 0;
       Max = [x_new y_new OptZ_new];
       Extreme = Max;
elseif abs(Hdet) > 0 && d2fdx2 < 0;
       Min = [x_new y_new OptZ_new];
       Extreme = Min;
elseif abs(Hdet) <0;
       Saddle = [x_new y_new OptZ_new];
       Extreme = Saddle;
end

%Plot for found mins and maxs
hold on
surf(X,Y,Z)
plot3(Extreme(1),Extreme(2),Extreme(3),'m*','MarkerSize',10)

Also, it seems to me that these two lines should be done once prior to the for loop (i.e., moved out of this section of code):

hold on
surf(X,Y,Z)

And also it seems that the rest of the code should be inside and at the end of the for loop, so you get four points plotted ... one for each iteration of the loop.

[Max]

I have been messing around with your suggestion and all that happens is the program writes it all to the max array which is nothing but negative values, lol, I do know the values I am looking for is z_max = 0.8484.

I am thinking that the program need to store the max/min/saddle values before the hessian but the hessian is what confirms the points as a max/min/saddle.

I'm going to try defining max = d2fdx2, those 2nd derivative values are already defined to include x_new and y_new.

[Max]

Still not working, OptZ_new is already defined as being equal to [x_new y_new] so I think that is just doubling it up

[Max]

The for loop finds the values and the hessian just confirms though, its subsequent to the values the for loop finds i think. I just need the hessian part to confirm them as a local min/max then store the point to plot.

[James Tursa]

OK, now you are asking for more than just plotting help. I haven't looked in detail at all of your code, but here are some obvious problems you need to fix:

You initialize these variables as 0's:

GS_Error = zeros(1,5); SA_Error = zeros(1,5);

And then test them as follows:

        while SA_Error>.01
                :
            while abs(GS_Error)>.01

First, testing a vector is not what you probably wanted ... I am guessing you wanted to use SA_Error(i) and GS_Error(i) here.

Second, these vectors start at 0's, so on the first pass they fail the test! You never get inside to execute the code you painstakingly wrote! SO maybe this instead:

GS_Error = ones(1,5); SA_Error = ones(1,5);

And I am not sure what to make of your indexing. Your for loop uses w but inside you use i. So something looks like it needs to be fixed there as well.

[Max]

it calls on w to iterate x_new and y_new through the array guess values defined as x_new = Guess(w,1), y_new = Guess(w,2)

I will mess around with what you have pointed out though, thank you

Answer 2

[Umair Asif]

Max were you able to complete the code? I am working on something similar to what you have displayed so by any chance you got that program to work? Plz respond asap

[Max]

Yes, I was able to get it to work. What are you stuck with?