%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%        MULTIVARIATE PLOTTING AND OPTIMIZATION               %%
%%        --------------------------------------               %%
%%   A brief analysis of a sample function f:R^2->R            %% 
%%   using contour and surface plotting, followed by           %% 
%%   a multivariate Newton algorithm for finding optimiziers.  %% 
%%   This may be easily modified to extend to functions of n   %% 
%%   variables.                                                %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% 1.  Create a grid and evaluate FF (:= f) for each point 
%%      on the grid, storing the results in the matrix A.
%% 
      x=-3:.2:3; y=x;
      n_x=size(x,2);n_y=size(y,2);
      for i=1:n_x;
        xx=x(1,i);
        for j=1:n_y;
         yy=y(1,j);
         A(j,i)=FF(xx,yy);
         end
        end
%%
%% 2. a) Compute the vectors to be used in the overlay quiver
%%       plot using the gradient function. 
%%
       [px,py]=gradient(A,.2,.2);
%%    b) Set up the vector v of contour values.
        v=-40:10:30; v=[v -8 -8.5 -6.5 -5.1852 -4.5];
%%    c) Contour, holding the plot and overlay the gradient field.
        CS=contour(A,v,x,y,'w'); hold on;
        pause(2);
        quiver(x,y,px,py);
%%    d) Create a new figure handle, and surface plot f over the 
%%       indicated domain.
        figure;
        mesh(x,y,A);
%% 3. Finally, isolate and obtain accurate approximations to 
%%    the critical points for f using Newton's method for  
%%    functions of 2 variables:
%%
        x_1=input('x_1 ='),x_2=input('x_2 =')
		while(x_1 ~= [] & x_2 ~= [])
          Newton_N(x_1,x_2)
          x_1=input('x_1 ='),x_2=input('x_2 =')
		  end
		
		
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%    Newton's Method for Optimizer of a Multivariate Function %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%  We use FF.m for the function, G.m for the gradient, and H.m
%%  for the Hessian. These need to be modified appropriately if 
%%  the objective function is to be modified.
  function [Newton_N] = Newton_N(x_1,x_2);
   A=[]; format long E
%%
%%   Select Parameters: initial guess & error tolerance
   x_current = [x_1,x_2]',  eps= 1.e-10
%%
%%   Initialize variables
%%   
   diff_current = (eps +1)*[1 1]';   % artificial to start
%%
%%   Begin interation loop
   while max(abs(diff_current)) > eps,
     grad  = G(x_current);          % Evaluate derivative  
     hess  = H(x_current);          % Evaluate 2-nd derivative 
     diff_new = -hess\grad;         % compute the step
     x_new = x_current + diff_new;  % update approx. minimizer
	 diff_current = diff_new;  
	 x_current = x_new ;       
     A=[ A                             % build matrix of results
          [ x_current'  diff_current' ]  ];     
	 end                               % end of do-while
   Newton_N= A;
   end
     
function [FF]=FF(x,y) 
%     x=z(1,1); y=z(2,1);
     FF= x*x*x + y*y*y -2*x*x +3*y*y-8;
	  end
function [G] = G(z) 
      x=z(1,1); y=z(2,1);
      G(1,1)=3*x*x-4*x; G(2,1)=3*y*y+6*y;
     end
function [H] = H(z) 
      x=z(1,1); y=z(2,1);
      H(1,1)=6*x-4; H(2,1)=0;
      H(1,2)=H(2,1); H(2,2)=6*y+6;
     end

x_current =[1 1]', eps = 1.000000000000000e-10
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                    Sample Output                        %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     x_current(1)       x_current(2)         diff(1)           diff(2)  
   ----------------   ----------------   ----------------  -----------------
   1.50000000000000   0.25000000000000   0.50000000000000  -0.75000000000000
   1.35000000000000   0.02500000000000  -0.15000000000000  -0.22500000000000
   1.33353658536585   0.00030487804878  -0.01646341463415  -0.02469512195122
   1.33333336430743   0.00000004646115  -0.00020322105842  -0.00030483158763
   1.33333333333333   0.00000000000000  -0.00000003097410  -0.00000004646115
   1.33333333333333   0.00000000000000  -0.00000000000000  -0.00000000000000