## Ejemplo 1: Encontrar los puntos de mayor y menor distancia al 
origen de la curva definida por x²+3xy+y²-4=0

kill(all)$

wxdraw2d(grid=true, line_type=solid, line_width=3, nticks=200,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
        color=blue, implicit(x^2+3*x*y+y^2-4=0, x, -3,3, y, -3,3), 
        line_width=1,  color=red, explicit(x,x,-3,3),  proportional_axes = xy )$

### Resolviendo a lo bestia

cond:solve(x^2+3*x*y+y^2-4=0,y);

define(fy1(x),rhs(cond[1]));

define(fy2(x),rhs(cond[2]));

define(f1(x),x^2+(fy1(x))^2);

define(f2(x),x^2+(fy2(x))^2);

/* Buscando minimos en la primera rama */
diff(f1(x),x);

der1:ratsimp(%);

solve(der1=0,x);

ls:num(der1)/3,ratsimp;

iz:first(ls);
de:ls-iz;

iz^2=de^2;

sol:solve(%,x);

sol[1];ev(der1,sol[1]);

sol[2];ev(der1,sol[2]);

ev(fy1(x),sol[1]);

diff(f1(x),x,2);

ev(%,sol[1]),ratsimp;

/* Buscando minimos en la segunda rama */
diff(f2(x),x);

der2:ratsimp(%);

ls:num(der2)/3,ratsimp;

iz:first(ls);
de:ls-%;

iz^2=de^2;

sol:solve(%,x);

ev(der2,sol[1]);

ev(der2,sol[2]);

ev(fy2(x),sol[2]);

diff(f2(x),x,2);

ev(%,sol[2]),ratsimp;

### Con los multiplicadores de Lagrange

kill(all)$

define(f(x,y),x^2+(y)^2);
define(L(x,y,l),f(x,y)+l*(x^2+y^2+3*x*y-4));
define(lig(x,y),x^2+y^2+3*x*y-4);

e1:diff(L(x,y,l),x);
e2:diff(L(x,y,l),y);
el:diff(L(x,y,l),l);

realonly:true$
sol:solve([e1,e2,el],[x,y,l]);
nn:length(%);

for kk:1 thru nn do print("f =",ev(f(x,y),sol[kk]) , "* Punto ", sol[kk] )$

define(hes(x,y,l),hessian(L(x,y,l),[x,y,l]))$
vec:matrix(  [dx],  [dy] )$ vect:[dx,dy]$

minor(hes(x,y,l),3,3)$
define(FQ0(dx,dy),expand(vect.%.vec));

k:1$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k]);
ratsimp(%);

k:2$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k]);
ratsimp(%);

ptos1:makelist( [rhs(sol[k][1]),rhs(sol[k][2])] , k,1,2);

wxdraw2d(grid=true, line_type=solid, line_width=3, nticks=200,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
        color=blue, implicit(lig(x,y)=0, x, -3,3, y, -3,3), 
        line_width=1,  color=red, explicit(x,x,-3,3),  proportional_axes = xy,
        point_type = filled_circle , point_size =2,   color=black, points(ptos1) )$

## Ejemplo 3: Encontrar la mínima distancia al origen  de la curva
  (x − 3)² + (y − 4)² = a² con a=10 y a=2

kill(all)$

a:10$
define(f(x,y),x^2+(y)^2);
define(lig(x,y),(x-3)^2+(y-4)^2-a^2);
define(L(x,y,l),f(x,y)+l*lig(x,y));

wxdraw2d(grid=true, line_type=solid, line_width=3, nticks=200,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
        color=blue, implicit(lig(x,y)=0, x, -10,15, y, -10,15), 
        proportional_axes = xy , 
     point_type = filled_circle , point_size =2,   color=black, points([[0,0]]) )$

realonly:true$
e1:diff(L(x,y,l),x);
e2:diff(L(x,y,l),y);
el:diff(L(x,y,l),l);
sol:solve([e1,e2,el],[x,y,l]);
nn:length(%);

for kk:1 thru nn do print("f =",ev(f(x,y),sol[kk]) , "* Punto ", sol[kk] )$

define(hes(x,y,l),hessian(L(x,y,l),[x,y,l]));
vec:matrix( [dx],  [dy] )$
vect:[dx,dy]$
minor(hes(x,y,l),3,3);
define(FQ0(dx,dy),expand(vect.%.vec));

k:1$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k]);

k:2$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k]);

ptosmax:makelist( [rhs(sol[k][1]),rhs(sol[k][2])] , k,1,1);
ptosmin:makelist( [rhs(sol[k][1]),rhs(sol[k][2])] , k,2,2);

wxdraw2d(grid=true, line_type=solid, line_width=2, nticks=200,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
        color=blue, implicit(lig(x,y)=0, x, -10,15, y, -10,15), 
        proportional_axes = xy , point_type = filled_circle , point_size =2,   
        color=black, points([[0,0]]), 
        color=red, points(ptosmax) , 
        color=dark_green, points( ptosmin))$

## Ejemplo 2. Encontrar los extremos de f(x,y) = xy con x²+y²=1

kill(all)$
define(f(x,y),x*y);
define(lig(x,y),x^2+y^2-1);
define(L(x,y,l),f(x,y)+l*lig(x,y));

realonly:true$
eq1:diff(L(x,y,l),x);
eq2:diff(L(x,y,l),y);
eq3:diff(L(x,y,l),l);
sol:solve([eq1,eq2,eq3],[x,y,l]);
nn:length(%);

eq1*y-eq2*x,ratsimp;

for kk:1 thru nn do print("f =",ev(f(x,y),sol[kk]) , "* Punto ", sol[kk] )$

define(hes(x,y,l),hessian(L(x,y,l),[x,y,l]));
vec:matrix( [dx],  [dy] )$
vect:[dx,dy]$
minor(hes(x,y,l),3,3);
define(FQ0(dx,dy),expand(vect.%.vec));

k:1$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k]);

k:2$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k]);

k:3$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k]);

k:4$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k]);

ptos1:makelist( [rhs(sol[k][1]),rhs(sol[k][2]),ev(f(x,y),sol[k])] , k,1,2);
ptos2:makelist( [rhs(sol[k][1]),rhs(sol[k][2]),ev(f(x,y),sol[k])] , k,3,4);

wxdraw3d( xu_grid = 25, yv_grid = 25, view = [78,15], surface_hide   = true,
             parametric_surface( cos(u), sin(u), sin(u)*cos(u), u, 0, 2*%pi, v, -%pi, %pi) ,
             point_type = filled_circle , point_size =2, 
             color=red, points(ptos1),  color=black, points(ptos2)  )$

## Ejemplo 4: Extremos de  x²+y²-x-y+1 en la región  x²+y²<=a, a=2 (a=1)

kill(all)$
define(f(x,y),x^2+y^2-x-y+1);

### Extremos libres

eq1:diff(f(x,y),x);
eq2:diff(f(x,y),y);

sol1:solve([eq1,eq2],[x,y]);

ev(f(x,y),sol1);

define(hes(x,y),hessian(f(x,y),[x,y]));

### Extremos en la frontera

define(lig(x,y),x^2+y^2-2);
define(L(x,y,l),f(x,y)+l*lig(x,y));

e1:diff(L(x,y,l),x);
e2:diff(L(x,y,l),y);
el:diff(L(x,y,l),l);
sol:solve([e1,e2,el],[x,y,l]);
nn:length(%);

for kk:1 thru nn do print("f =",ev(f(x,y),sol[kk]) , "* Punto ", sol[kk] )$

 print("f =",ev(f(x,y),sol1) , "* Punto ", sol1 )$

define(hes(x,y,l),hessian(L(x,y,l),[x,y,l]));
vec:matrix( [dx],  [dy] )$
vect:[dx,dy]$

minor(hes(x,y,l),3,3);
define(FQ0(dx,dy),expand(vect.%.vec));

k:1$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k]);
ratsimp(%);

k:2$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k]);
ratsimp(%);

wxdraw3d( xu_grid = 20, yv_grid = 20, view = [100,35],   
          point_type=filled_circle, point_size = 2,color = red, 
          points([[1,1,f(1,1)]]),
          surface_hide   = true, color = blue,
          parametric_surface( sqrt(2)*cos(u), sqrt(2)*sin(u),
          2*sin(u)^2+2*cos(u)^2-sqrt(2)*cos(u)- sqrt(2)*sin(u)+1, u, 0, 2*%pi, v, -%pi, %pi),
         color =dark_green,  points([[-1,-1,f(-1,-1)]])
    )$

wxdraw3d( xu_grid = 20, yv_grid = 20, view = [100,35],   
          point_type=filled_circle, point_size = 2,color = red, 
          points([[1,1,f(1,1)]]),
          surface_hide   = true, color = blue,
          parametric_surface( sqrt(2)*cos(u), sqrt(2)*sin(u),
          2*sin(u)^2+2*cos(u)^2-sqrt(2)*cos(u)- sqrt(2)*sin(u)+1, u, 0, 2*%pi, v, -%pi, %pi),
         color =dark_green,  points([[-1,-1,f(-1,-1)]]),  color =black,  points([[1/2,1/2,f(1/2,1/2)]])
    )$

¿Es realmente sol[2] un mínimo local?

sol[2]; ev(f(x,y),sol[2]);
k:2$
pmin:[rhs(sol[k][1]),rhs(sol[k][2])];

a:1$ b:1$
wxdraw2d(grid=true, line_type=solid, line_width=3, nticks=200,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
        color=blue, implicit(lig(x,y)=0, x, a-3,a+1, y,b-3,b+1), 
        proportional_axes = xy , color=black, explicit(x,x,0,1),
        point_type = filled_circle , point_size =2,   color=red, points([pmin]) )$

wxdraw2d(yrange=[0.65,1.1], xrange=[a-.2,a+.02],
    point_type = filled_circle,  point_size = 2, line_type=solid, line_width = 3,
    color = red, points([[a,f(a,a)]]), 
    color =dark_green, explicit(f(x,x),x,a-.2,a), /*  me acerco por y=x */
    color=blue, explicit(f(x,sqrt(2-x^2)),x,a-.2,a)  /*  me acerco por la frontera */  );

wxdraw3d( xu_grid = 20, yv_grid = 20, view = [83,57], surface_hide = true, dimensions = [750,750],  
    line_type=solid, line_width=3, xrange=[0.25,1.5], yrange=[0.25,1.5], zrange=[0.5,1.5], 
    parametric_surface( sqrt(2)*cos(u), sqrt(2)*sin(u), 2*sin(u)^2+2*cos(u)^2-sqrt(2)*cos(u)- 
        sqrt(2)*sin(u)+1, u, 0, %pi/2, v, 0, %pi) , point_type=filled_circle, point_size =2,
    color = red, points([[a,b,f(a,b)]]))$

wxdraw3d( xu_grid = 20, yv_grid = 20,view = [83,57], surface_hide = true,dimensions = [750,750],  
    line_type=solid, line_width=3, xrange=[0.25,1.5], yrange=[0.25,1.5],  zrange=[0.5,1.5], 
    parametric_surface( sqrt(2)*cos(u), sqrt(2)*sin(u), 2*sin(u)^2+
        2*cos(u)^2-sqrt(2)*cos(u)-sqrt(2)*sin(u)+1, u, 0, %pi/2, v, 0, %pi), 
    point_type=filled_circle, point_size =2,color = red, points([[a,b,f(a,b)]]),  
    color =dark_green,  parametric_surface( t, t,  f(t,t), t, 0, 1, t, 0,1 ) )$

Que ocurre en el caso que impongamos x>=0. Luego no tenemos máximo global (estaba en x<0).
Los calculos muestran 

f(0,y);

wxplot2d([f(0,y)],[y,-sqrt(2),sqrt(2)])$

Punto (0,1/2) no es mínimo local

f(0,1/2),float;

f(0.01,1/2);

f(0,1/2)-f(x,1/2),factor;

/* PUNTO (0,1/2) Me acerco por las x, no hay mínimo */ 
wxplot2d([f(x,1/2)],[x,0,.5]);

Punto (0,-sqrt(2)) 

f(0,-sqrt(2));float(%);

Punto (0,sqrt(2))

f(0,sqrt(2));float(%);

define(rec(x,t),+sqrt(2)-sqrt(2)*(x)/t);

wxdraw2d(grid=true, line_type=solid, line_width=3, nticks=200, 
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
        color=blue, implicit(lig(x,y)=0, x, a-3,a+1, y,b-3,b+1), 
     color=black, explicit(rec(x,1.3),x,0,1), explicit(rec(x,1),x,0,1),explicit(rec(x,.6),x,0,1),
        point_type = filled_circle , point_size =2, color=red, points([[0,sqrt(2)]]) )$

rec(0,0.6);

xx:0.001$
f(0,sqrt(2)),float;
f(xx,rec(xx,0.6)),float;
f(xx,rec(xx,1.3)),float;
f(xx,rec(xx,1.0)),float;

f(x,rec(x,a))-f(0,sqrt(2)),factor;

wxplot2d( f(x,rec(x,1))-f(0,sqrt(2))  ,[x,0,.1] );

## Ejemplo 5 (dos ligaduras) Hayar los puntos de mayor y menor distancia al origen de la curva definida por las ecuaciones x²+y²-z²-xy-1=0 y-1=0.

kill(all);

define(f(x,y,z),x^2+y^2+z^2);

define(lig1(x,y,z),x^2+y^2-z^2-x*y-1);
define(lig2(x,y,z),x^2+y^2-1);

jacobian([lig1(x,y,z),lig2(x,y,z)],[x,y,z]);

rank(%);

wxdraw3d(color=green,  view = [77,237], x_voxel = 15,
       y_voxel = 15, z_voxel = 15, enhanced3d = false, color= blue, 
implicit(x^2+y^2-1=0,x,-8,8,y,-8,8,z,-12,12),     
color = cyan, implicit(lig1(x,y,z)=0,x,-8,8,y,-8,8,z,-12,12),color = red,
point_type=filled_circle, point_size = 2,points([[0,0,0]])
)$

 wxdraw3d( xu_grid = 20, yv_grid = 20, view = [95,120], surface_hide = true,
             parametric_surface( cos(u), sin(u), sqrt(-sin(u)*cos(u)), u, %pi/2, %pi, v, %pi/2, %pi) ,
             parametric_surface( cos(u), sin(u), -sqrt(-sin(u)*cos(u)), u, %pi/2, %pi, v, %pi/2, %pi),
             parametric_surface( cos(u), sin(u), sqrt(-sin(u)*cos(u)), u, 3*%pi/2, 2*%pi, v, 3*%pi/2, 2*%pi) ,
             parametric_surface( cos(u), sin(u), -sqrt(-sin(u)*cos(u)), u, 3*%pi/2, 2*%pi, v, 3*%pi/2, 2*%pi), 
      point_type=filled_circle, point_size = 1.5, color = red,points([[0,0,0]]) );

define(L(x,y,z,l,m),f(x,y,z)+l*lig1(x,y,z)+m*lig2(x,y,z));

e1:diff(L(x,y,z,l,m),x);
e2:diff(L(x,y,z,l,m),y);
e3:diff(L(x,y,z,l,m),z);
el:diff(L(x,y,z,l,m),l);
em:diff(L(x,y,z,l,m),m);

realonly:true$
sol:solve([e1,e2,e3,el,em],[x,y,z,l,m])$
nn:length(%);

for kk:1 thru nn do print("f =",ev(f(x,y,z),sol[kk]) , "* Punto ", sol[kk] )$

/* Con el hessiano */;

define(hes(x,y,z,l,m),hessian(L(x,y,z,l,m),[x,y,z,l,m]));
vec:matrix( [dx],  [dy],  [dz] )$
vect:[dx,dy,dz]$
M0:minor(minor(hes(x,y,z,l,m),5,5),4,4);
define(FQ0(dx,dy,dz),expand(vect.M0.vec));

eql1:jacobian([lig1(x,y,z)],[x,y,z]).[dx,dy,dz];
eql2:jacobian([lig2(x,y,z)],[x,y,z]).[dx,dy,dz];
difligy:solve([eql1,eql2],[dx,dz])[1];
difligx:solve([eql1,eql2],[dy,dz])[1];

/* Distintas soluciones  */
k:1$
sol[k];

e1:ev(eql1,sol[k]);
e2:ev(eql2,sol[k]);
sollig:solve([e1,e2],[dx,dy])[1];

ev(FQ0(dx,dy,dz),[sollig,sol[k]]);

/* Distintas soluciones  */
k:7$
sol[k];

e1:ev(eql1,sol[k]);
e2:ev(eql2,sol[k]);
sollig:solve([e1,e2],[dx,dz])[1];

ev(FQ0(dx,dy,dz),[sollig,sol[k]]);

p1:makelist([rhs(sol[k][1]),rhs(sol[k][2]),rhs(sol[k][3]) ],k,1,4); /* min */
p2:makelist([rhs(sol[k][1]),rhs(sol[k][2]),rhs(sol[k][3]) ],k,5,8); /* max */

 wxdraw3d( xu_grid = 20, yv_grid = 20, view = [90,120], surface_hide = true,
             parametric_surface( cos(u), sin(u), sqrt(-sin(u)*cos(u)), u, %pi/2, %pi, v, %pi/2, %pi) ,
             parametric_surface( cos(u), sin(u), -sqrt(-sin(u)*cos(u)), u, %pi/2, %pi, v, %pi/2, %pi),
             parametric_surface( cos(u), sin(u), sqrt(-sin(u)*cos(u)), u, 3*%pi/2, 2*%pi, v, 3*%pi/2, 2*%pi) ,
             parametric_surface( cos(u), sin(u), -sqrt(-sin(u)*cos(u)), u, 3*%pi/2, 2*%pi, v, 3*%pi/2, 2*%pi), 
             point_type=filled_circle, point_size = 1.5, color = red,points([[0,0,0]]), 
            color = black, points(p1),  color = dark_green, points(p2)  );

## Ejemplo 6: Encontrar los todos los puntos críticos y los posibles extremos de  x²+y²+z²-2x+4y-3 en el dominio  x²+y²+z²<=9  (¿y si x>=0?)

### Extremos libres

kill(all)$
define(f(x,y,z),(x)^2+y^2+(z)^2-2*x+4*y-3);

e1:diff(f(x,y,z),x);
e2:diff(f(x,y,z),y);
e3:diff(f(x,y,z),z);

realonly:true$
solf:solve([e1,e2,e3],[x,y,z]);

ev(f(x,y,z),solf);

hessian(f(x,y,z),[x,y,z]);

### Extremos en la frontera x^2+y^2+z^2=9

define(lig(x,y,z),x^2+y^2+z^2-9);

define(L(x,y,z,l),f(x,y,z)+l*lig(x,y,z));

e1:diff(L(x,y,z,l),x);
e2:diff(L(x,y,z,l),y);
e3:diff(L(x,y,z,l),z);
el:diff(L(x,y,z,l),l);

sol:solve([e1,e2,e3,el],[x,y,z,l]);
nn:length(%);

for kk:1 thru nn do print("f =",ev(f(x,y,z),sol[kk]) , "* Punto ", sol[kk] )$

print("f =",ev(f(x,y,z),solf) , "* Punto ", solf )$

ev(f(x,y,z),sol[1]),float;

ev(f(x,y,z),sol[2]),float;

Weierstrass: f( sol[1] ) >  f( sol[2] ) -> sol[1] maximo sobre la frontera y sol[2] minimo sobre la frontera

Combrobando los extremos con la segunda diferencial

define(hes(x,y,z,l),hessian(L(x,y,z,l),[x,y,z,l]));
vec:matrix( [dx],  [dy],  [dz] )$
vect:[dx,dy,dz]$
M0:minor(hes(x,y,z,l),4,4);
define(FQ0(dx,dy,dz),expand(vect.M0.vec));

kk:1$
jacobian([lig(x,y,z)],[x,y,z]);
ev(%,sol[kk])$
%.vect; sollig:solve(%,dy);
ev(FQ0(dx,dy,dz),sollig);
ev(%,sol[kk]),factor;

kk:2$
jacobian([lig(x,y,z)],[x,y,z]);
ev(%,sol[kk])$
%.vect; sollig:solve(%,dy);
ev(FQ0(dx,dy,dz),sollig);
ev(%,sol[kk]),factor;

¿Es realmente sol[2] un mínimo local? Vemos que pasa si nos acercamos por el interior de la región. 

sol[2];
float(sol[2]);
fmin:ev(f(x,y,z),%);

x0:rhs(sol[2][1])$ y0:rhs(sol[2][2])$

Nos acercamos a sol[2] por el plano z=0. La ligadura es entonces x^2+y^2=9

f(x,y,0);

lig(x,y,0);

wxdraw2d(  line_type=solid, line_width=3,   color =blue,  xlabel="x", ylabel="f(x,y,0) por la frontera",
                explicit(f(x,-sqrt(9-x^2),0),x,x0-.2,x0+.2) ,  
                  point_type=filled_circle, point_size=2,  color = dark_green, 
                  points([[x0,fmin]]) )$

/* Me acerco por las y con la x=x0 fija */
wxdraw2d(  line_type=solid, line_width=3,    xlabel="y", ylabel="f(x0,y,0), x fija",
             color =blue,  explicit( f(x0,y,0) ,y,y0-.2,y0) ,  
                  point_type=filled_circle, point_size=2,  color = dark_green, 
                  points([ [y0,f(x0,y0,0)]]) )$

/* Me acerco por las x con la y=y0 fija */
wxdraw2d(  line_type=solid, line_width=3,   xlabel="x", ylabel="f(x,y0,0), y fija",
               color =blue,  explicit( f(x,y0,0) ,x,x0-.2,x0) ,  
                  point_type=filled_circle, point_size=2,  color = dark_green, 
                points([  [x0,f(x0,y0,0)] ]   ) )$

### Extremos en la región x^2+y^2+z^2<=9 y x>=0

Imponiendo ahora que x>=0. Lo anterior se mantiene pero no tenemos máximo global pues este estaba en x<0. 
Luego el máximo ha de estar en la cara x=0. Estudiamos el problema de extremos si x=0.

kill(g)$
define(g(y,z),f(0,y,z));

/* Extremos libres */
eq1:diff(g(y,z),z);
eq2:diff(g(y,z),y);

sol1:solve([eq1,eq2],[y,z]);

yy0:rhs(sol1[1][1]), float; zz0:rhs(sol1[1][2]), float; 

g(y,z);

g(yy0,zz0);

define(hes(y,z),hessian(g(y,z),[y,z]));

¿Es un mínimo local de f?

print("f =",f(0,yy0,zz0) , "* Punto ", [0,yy0,zz0] )$

/* Restringido a la frontera x=0 */
wxdraw3d( xu_grid = 25, yv_grid = 25, view = [78,15],
    explicit(g(y,z),y,yy0-.1,yy0+.1,z,zz0-.1,zz0+.1),
    point_type=filled_circle, point_size = 2,color = red, 
    points([[yy0,zz0,g(yy0,zz0)]]))$

/* nos acercamos por las "y" y es mínimo */  
wxdraw2d(  line_type=solid, line_width=3,   xlabel="y", ylabel="f(0,y0,0) y=y0",
               color =blue,  explicit( f(0,y,0) ,y,yy0-.5,yy0+.5),  
                  point_type=filled_circle, point_size=2,  color = dark_green, 
                points([  [yy0,f(0,yy0,0)] ]   ), 
    color=black,label([concat("Punto ", "[", 0,",",yy0,",",zz0,"]",  ", f =", f(0,yy0,zz0) ) ,-2,-6.775] ) ) $

/* nos acercamos por las z  y es mínimo*/
wxdraw2d(  line_type=solid, line_width=3,   xlabel="z", ylabel="f(0,yy0,z) y=y0",
               color =blue,  explicit( f(0,yy0,z) ,z,zz0-.5,zz0+.5),  
               point_type=filled_circle, point_size=2,  color = dark_green, 
               points([ [zz0,f(0,yy0,zz0)] ]),
     color=black,label([concat("Punto ", "[", 0,",",yy0,",",zz0,"]",  ", f =", f(0,yy0,zz0) ) , 0,-6.775] ) )$

/* nos acercamos por el interior las x>0, no es mínimo */
xx0:0$
wxdraw2d(  line_type=solid, line_width=3,   xlabel="x", ylabel="f(x,yy0,0) x=x0",
               color =blue,  explicit( f(x,yy0,zz0) ,x,xx0,xx0+.5),  
               point_type=filled_circle, point_size=2,  color = dark_green, 
               points([ [zz0,f(0,yy0,zz0)] ]),
     color=black,label([concat("Punto ", "[", 0,",",yy0,",",zz0,"]",  ", f =", f(0,yy0,zz0) ) , 0.3,-7.1] ) )$

Con x=0 y la ligadura de la frontera 

lig(x,y,z);

remfunction(L)$
define(L(y,z,l),g(y,z)+l*(lig(0,y,z)));

e1:diff(L(y,z,l),y);
e2:diff(L(y,z,l),z);
el:diff(L(y,z,l),l);

solg:solve([e1,e2,el],[y,z,l]);
nn:length(%);

for kk:1 thru nn do print("f =",ev(f(0,y,z),solg[kk]) , "* Punto ", sol[kk] )$

sol[1] maximo absoluto en la frontera y sol[2] mínimo absoluto en la frontera. 
Por Weiestrass sol[1] ha de ser el máximo absoluto de f en toda la región. 

¿Es un mínimoo local   un mínimo local?   

xx0:0$ yy0:ev(y,solg[2])$  zz0:ev(z,solg[2])$

/* interior de la frontera x=0*/
wxdraw2d(  line_type=solid, line_width=3,   xlabel="x", ylabel="f(x,yy0,0) x=x0",
               color =blue,  explicit( f(0,-sqrt(9-z^2),z) ,z,zz0-.5,zz0+.5),  
               point_type=filled_circle, point_size=2,  color = dark_green, 
               points([ [zz0,f(0,yy0,zz0)] ]),
     color=black,label([concat("Punto ", "[", 0,",",yy0,",",zz0,"]",  ", f =", f(0,yy0,zz0) ) , 0.,-5.85] ) )$

/* Interior de la región x>0, "y" y "z" fijas */
wxdraw2d(  line_type=solid, line_width=3,   xlabel="x", ylabel="f(x,yy0,zz0)",
               color =blue,  explicit( f(x,yy0,zz0) ,x,xx0,xx0+.5),  
               point_type=filled_circle, point_size=2,  color = dark_green, 
               points([ [zz0,f(0,yy0,zz0)] ]),
     color=black,label([concat("Punto ", "[", 0,",",yy0,",",zz0,"]",  ", f =", f(0,yy0,zz0) ) , 0.3,-6.1] ) )$

## Ejemplo 7 (dos ligaduras): Sea la curva definida por el paraboloide z= x²+y² y el plano x+y+z=12.  

### Encontrar el punto más alto y el más bajo.

kill(all)$
define(f(x,y,z),z);

define(lig1(x,y,z),x+y+z-12);
define(lig2(x,y,z),x^2+y^2-z);

solve(lig1(x,y,z),z);

ev(lig2(x,y,z),%);

jacobian([lig1(x,y,z),lig2(x,y,z)],[x,y,z]);

rank(%);

wxdraw3d(color=green,  view = [80,22],    x_voxel = 15,
       y_voxel = 15, z_voxel = 15, enhanced3d = false, color= cyan, 
implicit(x^2+y^2-z=0,x,-8,8,y,-8,8,z,-2,20),     color= blue, 
implicit(z+y+x−12=0,x,-8,8,y,-8,8,z,-2,20),surface_hide = true
)$

/*  En paramétricas */
x(u):=-1/2+5/sqrt(2)*cos(u); y(u):=-1/2+5/sqrt(2)*sin(u);
wxdraw3d( xu_grid = 20, yv_grid = 20, view = [120,35],   
          parametric_surface( x(u),y(u), -y(u)-x(u)+12, u, 0, 2*%pi, v, -%pi, %pi) )$

define(L(x,y,z,l,m),f(x,y,z)+l*lig1(x,y,z)+m*lig2(x,y,z));

e1:diff(L(x,y,z,l,m),x);
e2:diff(L(x,y,z,l,m),y);
e3:diff(L(x,y,z,l,m),z);
el:diff(L(x,y,z,l,m),l);
em:diff(L(x,y,z,l,m),m);

sol:solve([e1,e2,e3,el,em],[x,y,z,l,m]);
nn:length(%);

for kk:1 thru nn do print("f =",ev(f(x,y,z),sol[kk]) , "* Punto ", sol[kk] )$

/* Con el hessiano */;

define(hes(x,y,z,l,m),hessian(L(x,y,z,l,m),[x,y,z,l,m]));
vec:matrix( [dx],  [dy],  [dz] )$
vect:[dx,dy,dz]$
M0:minor(minor(hes(x,y,z,l,m),5,5),4,4);
define(FQ0(dx,dy,dz),expand(vect.M0.vec));

/* Distintas soluciones  */
k:1$
sol[k];

eql1:jacobian([lig1(x,y,z)],[x,y,z]).[dx,dy,dz];
eql2:jacobian([lig2(x,y,z)],[x,y,z]).[dx,dy,dz];
e1:ev(eql1,sol[k]);
e2:ev(eql2,sol[k]);
sollig:solve([e1,e2],[dx,dz])[1];

ev(FQ0(dx,dy,dz),[sollig,sol[k]]);

/* Distintas soluciones  */
k:2$
sol[k];

eql1:jacobian([lig1(x,y,z)],[x,y,z]).[dx,dy,dz];
eql2:jacobian([lig2(x,y,z)],[x,y,z]).[dx,dy,dz];
e1:ev(eql1,sol[k]);
e2:ev(eql2,sol[k]);
sollig:solve([e1,e2],[dx,dz])[1];

ev(FQ0(dx,dy,dz),[sollig,sol[k]]);

p1:makelist([rhs(sol[k][1]),rhs(sol[k][2]),rhs(sol[k][3]) ],k,1,1);
p2:makelist([rhs(sol[k][1]),rhs(sol[k][2]),rhs(sol[k][3]) ],k,2,2);

wxdraw3d( xu_grid = 20, yv_grid = 20, view = [120,35],   
  parametric_surface( x(u),y(u), -y(u)-x(u)+12, u, 0, 2*%pi, v, -%pi, %pi) ,
  point_type=filled_circle, point_size = 1.5,  
  color = black, points(p1),
  color = dark_green, points(p2))$

### Hallar los puntos de mayor y menor distancia al origen

remfunction(f)$
define(f(x,y,z),x^2+y^2+z^2);

define(lig1(x,y,z),x+y+z-12);
define(lig2(x,y,z),x^2+y^2-z);

jacobian([lig1(x,y,z),lig2(x,y,z)],[x,y,z]);

rank(%);

x(u):=-1/2+5/sqrt(2)*cos(u);
y(u):=-1/2+5/sqrt(2)*sin(u);
wxdraw3d( xu_grid = 20, yv_grid = 20, view = [120,35],   
          parametric_surface( x(u),y(u), -y(u)-x(u)+12, u, 0, 2*%pi, v, -%pi, %pi),
     point_type=filled_circle, point_size = 1.5,   color = red, points([[0,0,0]]) 
)$

define(L(x,y,z,l,m),f(x,y,z)+l*lig1(x,y,z)+m*lig2(x,y,z));

e1:diff(L(x,y,z,l,m),x);
e2:diff(L(x,y,z,l,m),y);
e3:diff(L(x,y,z,l,m),z);
el:diff(L(x,y,z,l,m),l);
em:diff(L(x,y,z,l,m),m);

realonly:true$
sol:solve([e1,e2,e3,el,em],[x,y,z,l,m]);
nn:length(%);

for kk:1 thru nn do print("f =",ev(f(x,y,z),sol[kk]) , "* Punto ", sol[kk] )$

/* Con el hessiano */;

define(hes(x,y,z,l,m),hessian(L(x,y,z,l,m),[x,y,z,l,m]))$
vec:matrix( [dx],  [dy],  [dz] )$
vect:[dx,dy,dz]$
M0:minor(minor(hes(x,y,z,l,m),5,5),4,4);
define(FQ0(dx,dy,dz),expand(vect.M0.vec));

/* Distintas soluciones  */
k:1$
sol[k];

eql1:jacobian([lig1(x,y,z)],[x,y,z]).[dx,dy,dz];
eql2:jacobian([lig2(x,y,z)],[x,y,z]).[dx,dy,dz];
e1:ev(eql1,sol[k]);
e2:ev(eql2,sol[k]);
sollig:solve([e1,e2],[dx,dz])[1];

ev(FQ0(dx,dy,dz),[sollig,sol[k]]);

/* Distintas soluciones  */
k:2$
sol[k];

eql1:jacobian([lig1(x,y,z)],[x,y,z]).[dx,dy,dz];
eql2:jacobian([lig2(x,y,z)],[x,y,z]).[dx,dy,dz];
e1:ev(eql1,sol[k]);
e2:ev(eql2,sol[k]);
sollig:solve([e1,e2],[dx,dz])[1];

ev(FQ0(dx,dy,dz),[sollig,sol[k]]);

p1:makelist([rhs(sol[k][1]),rhs(sol[k][2]),rhs(sol[k][3]) ],k,1,1);
p2:makelist([rhs(sol[k][1]),rhs(sol[k][2]),rhs(sol[k][3]) ],k,2,2);

wxdraw3d( xu_grid = 20, yv_grid = 20, view = [80,60],   
  parametric_surface( x(u),y(u), -y(u)-x(u)+12, u, 0, 2*%pi, v, -%pi, %pi) ,
  point_type=filled_circle, point_size = 1.5,   color = red, points([[0,0,0]]),
  color = black, points(p1),
  color = dark_green, points(p2))$

## Ejemplo 8: Extremos de f(x,y)=log(x²+y²) en (x-1)²+y²<= 4.

kill(all);

define(f(x,y),log(x^2+(y)^2));

e1:diff(f(x,y),x);
e2:diff(f(x,y),y);
sol:solve([e1,e2],[x,y]); /* MAL */

/* En la circunferencia */
x(u):=1+2*cos(u); y(u):=2*sin(u);
wxdraw3d( xu_grid = 20, yv_grid = 20, view = [100,35],   
          point_type=filled_circle, point_size = 2, 
          surface_hide   = true, color = blue,
          parametric_surface( x(u), y(u),  log( x(u)^2+y(u)^2 ), u, 0, 2*%pi, v, -%pi, %pi ));

define(lig(x,y),(x-1)^2+y^2-4);
define(L(x,y,l),f(x,y)+l*lig(x,y));

e1:diff(L(x,y,l),x);
e2:diff(L(x,y,l),y);
el:diff(L(x,y,l),l);
sol:solve([e1,e2,el],[x,y,l]);
nn:length(%);

e1*y-e2*x,factor;

for kk:1 thru nn do print("f =",ev(f(x,y),sol[kk]) , "* Punto ", sol[kk] )$

define(hes(x,y,l),hessian(L(x,y,l),[x,y,l]))$
vec:matrix( [dx],  [dy] )$
vect:[dx,dy]$

minor(hes(x,y,l),3,3)$
define(FQ0(dx,dy),expand(vect.%.vec));

k:1$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k])$
ratsimp(%);

k:2$
sol[k];
jacobian([lig(x,y)],[x,y]);
ev(%,sol[k]).vect;
dlig:solve(%,dx);

ev(FQ0(dx,dy),dlig);
ev(%,sol[k])$
ratsimp(%);

p1:makelist([rhs(sol[k][1]),rhs(sol[k][2]), f(rhs(sol[k][1]),rhs(sol[k][2])) ],k,1,1);
p2:makelist([rhs(sol[k][1]),rhs(sol[k][2]), f(rhs(sol[k][1]),rhs(sol[k][2])) ],k,2,2);

wxdraw3d( xu_grid = 20, yv_grid = 20, view = [100,5],   
          point_type=filled_circle, point_size = 2, surface_hide   = true, color = blue,
          parametric_surface( x(u), y(u),  log( x(u)^2+y(u)^2 ), u, 0, 2*%pi, v, -%pi, %pi ),  
    color = red, points(p1),   color = dark_green, points(p2) )$

/* x=1 */;

f(1,y);

wxdraw2d( color = brown, line_width  = 3, 
                  explicit(f(1,y),y,-2,2)  )$

f(1,2);
float(%);

f(x(u),y(u))$
trigsimp(%);

ev(f(x(u),y(u)),u=%pi/2);float(%);

ev(f(x(u),y(u)),u=%pi/2-.01);

wxdraw3d( xu_grid = 20, yv_grid = 20,view = [100,50], surface_hide = true, 
    line_type=solid, line_width=3,  parametric_surface( x(u), y(u),  log( x(u)^2+y(u)^2 ), u, -%pi/2,%pi/2, v, 0, %pi ),  
    point_type=filled_circle, point_size =2,  color = dark_green, points(p2) ,   color =black , points([[1,0,f(1,0)]]) , 
    color =red, points([[1,2,f(1,2)]]) ,points([[1,-2,f(1,-2)]]) ,
    color =red,  parametric_surface( 1, t,  f(1,t), t, -2, 2, t, -2,2 ) )$