kill(all)$

simple optimization f(x,y) : 'fixed step,fx,fy', part 1a ,version 1.0a
peter.vlasschaert@gmail.com, 21/03/2019

info :  (graphics)  'see below'
==============
http://people.kzoo.edu/barth/math214/SurfacePlaneLine.html
rem :
info :  (maximum)
===============
simple optimization f(x,y) : 'fixed step,fx,fy,hessian', part 1b ,version 1.0a.
result : maximun

ref: titele :optimization of chemical processes .
        auth : D.M.Himmelblau.

function : f(x,y)  = ?
================

p1:f(x,y):= 4 -x^2-y^2;

∂f(x,y)/∂x = ?
============

p1a:'diff(f(x,y),x) =diff(f(x,y),x);

p1b:'diff(f(x,y),y) =diff(f(x,y),y);

gradient  :  grad(f)  = ∂f(x,y)/∂x*ex + ∂f(x,y)/∂y*ey
=======

p1c:[rhs(p1a),rhs(p1b)];

info : ex =>  unit vector x-direction, component = ∂f(x,y)/∂x.
           ey =>  unit vector y-direction, component = ∂f(x,y)/∂y.

 plot  : f(x,y) = 4-x^2-y^2
=====================

show more : contour levels ( with this option).
-------------------------------------------------------------------------

set_plot_option ([gnuplot_preamble, "set cntrparam levels 20"])$

a) contour plot (level surfaces = cross section of a plane with z axis)
=========================================================
=>  2D  from 3D

p2a:wxcontour_plot (f(x,y), [x, -4, 4], [y, -4, 4])$

b) plot3d
==========
=> 3D from 3D

p2b:wxplot3d([f(x,y), 4, [x, -4, 4], [y, -4, 4]])$

legend : z = 4
                 z = f(x,y) =4 -x^2-y^2

rem : will be maxinimum for function ' wxplot3d'
' more mathematics statement or needed to show this.'

start point : search the maximum (x0,y0)
------------------------------------------------------------

p3a : x0 = 0.5;
p3b : y0 = -1.0;

f(x0,y0) = 
======

p4:'f(x0,y0) = f(rhs(p3a),rhs(p3b));

step value
=========

p4a: h = 0.12;

p5a:x=rhs(p3a) + rhs(p4a)*subst([ x=rhs(p3a)],rhs(p1a));
p5b:y=rhs(p3b) + rhs(p4a)*subst([ y=rhs(p3b)],rhs(p1b));

p6a:x=rhs(p5a) + rhs(p4a)*subst([ x=rhs(p5a)],rhs(p1a));
p6b:y=rhs(p5b) + rhs(p4a)*subst([ y=rhs(p5b)],rhs(p1b));

p7a:x=rhs(p6a) + rhs(p4a)*subst([ x=rhs(p6a)],rhs(p1a));
p7b:y=rhs(p6b) + rhs(p4a)*subst([ y=rhs(p6b)],rhs(p1b));

p7a:x=rhs(p6a) + rhs(p4a)*subst([ x=rhs(p6a)],rhs(p1a));
p7b:y=rhs(p6b) + rhs(p4a)*subst([ y=rhs(p6b)],rhs(p1b));

p8a:x=rhs(p7a) + rhs(p4a)*subst([ x=rhs(p7a)],rhs(p1a));
p8b:y=rhs(p7b) + rhs(p4a)*subst([ y=rhs(p7b)],rhs(p1b));

p9a:x=rhs(p8a) + rhs(p4a)*subst([ x=rhs(p8a)],rhs(p1a));
p9b:y=rhs(p8b) + rhs(p4a)*subst([ y=rhs(p8b)],rhs(p1b));

p10a:x=rhs(p9a) + rhs(p4a)*subst([ x=rhs(p9a)],rhs(p1a));
p10b:y=rhs(p9b) + rhs(p4a)*subst([ y=rhs(p9b)],rhs(p1b));

p11a:x=rhs(p10a) + rhs(p4a)*subst([ x=rhs(p10a)],rhs(p1a));
p11b:y=rhs(p10b) + rhs(p4a)*subst([ y=rhs(p10b)],rhs(p1b));

p12a:x=rhs(p11a) + rhs(p4a)*subst([ x=rhs(p11a)],rhs(p1a));
p12b:y=rhs(p11b) + rhs(p4a)*subst([ y=rhs(p11b)],rhs(p1b));

p13a:x=rhs(p12a) + rhs(p4a)*subst([ x=rhs(p12a)],rhs(p1a));
p13b:y=rhs(p12b) + rhs(p4a)*subst([ y=rhs(p12b)],rhs(p1b));

p14a:x=rhs(p13a) + rhs(p4a)*subst([ x=rhs(p13a)],rhs(p1a));
p14b:y=rhs(p13b) + rhs(p4a)*subst([ y=rhs(p13b)],rhs(p1b));

p15a:x=rhs(p14a) + rhs(p4a)*subst([ x=rhs(p14a)],rhs(p1a));
p15b:y=rhs(p14b) + rhs(p4a)*subst([ y=rhs(p14b)],rhs(p1b));

p16a:x=rhs(p15a) + rhs(p4a)*subst([ x=rhs(p15a)],rhs(p1a));
p16b:y=rhs(p15b) + rhs(p4a)*subst([ y=rhs(p15b)],rhs(p1b));

p17a:x=rhs(p16a) + rhs(p4a)*subst([ x=rhs(p16a)],rhs(p1a));
p17b:y=rhs(p16b) + rhs(p4a)*subst([ y=rhs(p16b)],rhs(p1b));

p18a:x=rhs(p17a) + rhs(p4a)*subst([ x=rhs(p17a)],rhs(p1a));
p18b:y=rhs(p17b) + rhs(p4a)*subst([ y=rhs(p17b)],rhs(p1b));

go further  = >  limit will be  (0,0)

limit calculation
=============

pa:'limit(f(x,y), x, 0) = limit(f(x,y), x, 0) ;

pb:'limit(rhs(pa), y, 0) =limit(rhs(pa), y, 0);

Maximum
========

pc:subst([x=0,y=0],f(x,y));

code : max or min.
---------------------------
==========================
hh                 = step value
start_value = x0,y0  for search 
ff                    = f(x,y)
nn                  = number of steps
==========================

fxy_max_or_min(hh,xx0,yy0,ff,nn):= block(
[], h: hh, f1:ff,x0:xx0,y0:yy0,n:nn,
    fx:diff(f1,x), fxx:subst([x=x0,y=y0],fx),
    fy:diff(f1,y), fyy:subst([x=x0,y=y0],fy),    
            xx : x0+h*fxx,
            yy : y0 +h*fyy,    
      for i:1  thru n
             do(
    xx:xx+subst([x=xx,y=yy],fx)*h,
    yy:yy+subst([x=xx,y=yy],fy)*h
                         ),
    f11:subst([x=xx,y=yy],f1),
           [xx,yy,f11]   
)$

example : 'code'
========

pd:fxy_max_or_min(0.12,0.5,-1.0,4-x^2-y^2,10);

list : output:  'pd'
==============

x-value : maximum
-----------------------------

pd1:pd[1];

y-value : maximum
-----------------------------

pd2:pd[2];

z-value : maximum
-----------------------------

pd3:pd[3];

'-------------------------------------------------------------------------------------------------------------------------------'

rem : n = 10  => same result  'p15a,p15b'
rem : wx , inside document for plot3d, contour_plot.

'--------------------------------------------------------------------------------------------------------------------------------'