# Ejemplos de aplicaciones del cálculo diferencial
a funciones de varias variables

## La diferencial en un campo vectorial

Utilizamos esa denominación para las funciones de R^n (o un subconjunto abierto de R^n) en R^m

En cuanto al comando diff(F(x,y)) da igual que F sea una función escalar o vectorial

    print("las funciones coordenadas y sus diferenciales")$
f1(x,y):=x*y;
diff(f1(x,y));
f2(x,y):= x^2+y^3;
diff(f2(x,y));
f3(x,y):=cos(x*y);
diff(f3(x,y));

    print("la diferencial del campo vectorial: todas de golpe")$
f(x,y):= [f1(x,y),f2(x,y),f3(x,y)];
expand(f(x,y));
diff(%);

## Plano tangente y gradiente

    a:1/3$ b:1/3$                     /* punto seleccionado */
f(x,y):= -(x^2+y^2);                  /* superficie         */
jacobian([f(x,y)],[x,y]);             /* matriz de las derivadas parciales */
gradiente:ev(%, [x=a,y=b]);           /* valor en (a,b)  */
g(x,y):=f(a,b)+gradiente.[x-a,y-b];   /* plano tangente en (a,b) */



wxdraw3d( surface_hide=true, view=[75,25],
    explicit(f(x,y),x,-2,2,y,-2,2),
    color=red,
    explicit( g(x,y), x, -3/2,3/2, y, -3/2,3/2))$

## Diferenciabilidad, máximos y mínimos

        ENUNCIADO DEL EJERCICIO 1

Los polinomios y el seno son funciones continuas. En consecuencia, salvo en el origen,
que requiere un estudio más específico, f es una función continua. 
Si la dibujamos, parece que la continuidad en (0,0) está garantizada.  

f(x,y):= sin(x^2+y^2)/(x^2+y^2);
wxdraw3d( enhanced3d=true, 
          colorbox=false,
          palette=[6,3,6],
    explicit ( f(x,y), x,-3.5,3.5, y,-3.5,3.5))$

La continuidad en el origen requiere probar que lim f(x,y) =1 
cuando (x,y) tiende a (0,0). 

/* Por sí mismo, Maxima sólo es capaz de calcular límites en una variable */

? limit;

/* En este caso, con eso basta, porque haciendo un cambio de variable
a coordenadas polares, que (x,y) tienda a 0 equivale a que r tiende a 0  */

subst( x^2+y^2=r^2, f(x,y));

limit(%,r,0);

print("Esto prueba la continuidad de f en (0,0), que muestra la gráfica")$

Al igual que ocurre con la continuidad, la diferenciabilidad, salvo en (0,0), es consecuencuencia
de teoremas generales: tanto los polinomios como el seno son funciones infinitamente derivables
y en consecuencia lo mismo ocurre con f. 

print("Salvo en el origen, que de momento excluimos, la diferencial está dada por")$                 
diff(f(x,y));    /*  del(x) y del(y)  representan a  dx  y dy en la notación clásica  */   
fullratsimp(%); 

print("que podemos expresarla también en términos del gradiente")$

jacobian([f(x,y)],[x,y]);
gradiente: fullratsimp(%);
print("Las derivadas parciales son funciones continuas. 
Por tanto, f es de clase C^1(R^2) y en consecuencia es diferenciable")$

Pasamos ahora a analizar la diferenciabilidad en el (0,0). 
Comenzamos estudiando la existencia de la derivada en la dirección del vector unitario genérico v: [a,b]

(f(0+t*a,0+t*b)-1)/t;
limit(%,t,0);

print("En consecuencia, todas las derivadas direccionales salen 0.
La diferencial, caso de existir también ha de ser 0.
La imagen parece indicar que en (0,0) hay plano tangente a la superficie que es horizontal.
Esa es justamente la traducción geométrica de lo anterior.")$

print("Que la diferencial de f en el origen exista y valga 0 
equivale a que limit( (f(h,k)-1)/sqrt(h^2+k^2)=0 cuando (h,k) tiende a cero.
Procediendo como más arriba hacemos un cambio de variable a polares")$
(f(h,k)-1)/sqrt(h^2+k^2);
subst(h^2+k^2=r^2, %);
limit(%,r,0);
print("En resumen, f es continua y diferenciable en todos los puntos")$

El cálculo de los extremos absolutos es muy sencillo ya que, escrita en coordenadas polares,
f se expresa como sin(r^2)/r^2 y podemos dibujar su gráfica

wxdraw2d( explicit(sin(r^2)/r^2, r, 0, 10))$
print("La gráfica muestra que la función tiene un máximo absoluto para r=0 
y un mínimo absoluto cerca de r=2, que vamos a tratar de calcular");

ec: diff( sin(r^2)/r^2, r);
solve(%=0,r);
print("Como no puede resolver la ecuación de forma exacta, 
le vamos a ayudar a encontrar una solución aproximada, inspirándonos en la gráfica")$

find_root(ec=0, r,2,2.5);
ev(sin(r^2)/r^2, r=%);

             ENUNCIADO DEL EJERCICIO 2

   /* DEFINIMOS g Y LA DISTANCIA A P DE UN PUNTO CUALQUIERA */
g(x,y):=x^2+y^2;         
D(x,y,z):=(x-1)^2+y^2+(z-2)^2;

           CONDICIÓN NECESARIA DE EXTREMO

Los puntos (x,y,z) pertenecientes a la gráfica de g cumplen que z=g(x,y)

G(x,y):=subst(z=g(x,y),D(x,y,z)); /* distacia de P a un punto de g */
G(x,y);                           /* ver el resultado de la sustitución */
print("G es la función a minimizar. 
Desarrollamos y calculamos su diferencial (=) jacobiana")$
expand(G(x,y));
J: jacobian([G(x,y)],[x,y]);


print("Calculamos los puntos donde la diferencial (J)  se anula,
lo cual significa que sean nulas las derivadas parciales respecto: x e y ")$
J[1];                   /* no proporciona la primera coordenada es J[1,1]  */
J[1,1];                 /* es la derivada parcial respecto a x */
solve( [J[1,1],J[1,2]],  [x,y]); /* no es necesario poner J(1,*)=0  */;

print("Maxima ha encontrado tres puntos posibles de extremo;
podemos comprobarlo por otro procedimiento,
le decimos que factorice y resolvemos")$
factor(J);

Claramente y=0 es una solución de la segunda condición y con ella 2*x^3-3*x-1=0 
lo que implica resolver la ecuación siguiente

solve(2*x^3-3*x-1,x);

La otra posibilidad es 2*x^3-3*x-1=0, pero entonces la primera coordenada es -2.
Moraleja, Maxima no se ha equivocado en la determinación de los puntos críticos.

          CONDICIÓN SUFICIENTE DE EXTREMO 

Para analizar si los puntos críticos corresponden, o no, a extremos 
hay que analizar la diferencial segunda, es decir la hessiana, en esos puntos.

H:hessian([G(x,y)],[x,y]);   /* definimos H */
H;                           /* le pedimos que nos muestre su valor */

print("Procedemos a evaluarla en cada uno de los puntos críticos")$

Q:ev(H,x=(1-sqrt(3))/2, y=0);
%,numer;

print("Forma definida negativa: máximo relativo estricto en
    x=(1-sqrt(3))/2, y=0 ")$

Q:ev(H,x= (sqrt(3)+1)/2, y=0);
%,numer;
print("Forma definida positiva: mínimo relativo estricto en
    x=(1+sqrt(3))/2, y=0 ")$

print("Calculamos esa distancia mínima")$
ev(G(x,y),  x=(1-sqrt(3))/2, y=0),expand;

Q:ev(H,x= -1, y=0);
print("Forma indefinida: punto de silla en
    x=-1, y=0 ")$

wxdraw3d(view=[50,320], key_pos=bottom_right,
    xtics=false, ytics=false, ztics=false, xlabel="",
    explicit(g(x,y),x,-2.5,2.5,y,-0.5,1.5),
color=black,point_size=2,point_type=filled_circle, key="P",
points([[1,0,2]]),
color=green, key="mínima distancia",
points([[(1+sqrt(3))/2,0,g((1+sqrt(3))/2,0)]]),
color=red, key="máxima distancia",
points([[(1-sqrt(3))/2,0,g((1-sqrt(3))/2,0)]]),
color=magenta, key="silla", point_type=circle,
points([[-1,0,g(-1,0)]])    
)$

## Extremos condicionados

    ENUNCIADO DEL EJERCICIO 3

f(x,y):=(x-1/8)^2+(y+4/5)^2;    /* función a optimizar */
g1(x,y):=x^2+(y+1/4)^2-1;                /* ligadura para item 1 */
g2(x,y):=y^4+2*x^2*y^2+2*y^2+x^4-2*x^2;  /* ligadura para item 2 */;

Visualizamos los correspondientes gráficos para ambas ligaduras
y el punto de referencia 

wxdraw2d( line_width=3, nticks=10000, proportional_axes=xy,
    key="g1",
implicit(g1(x,y), x,-2,2,y,-2,2),
        color=red, key="g2",
implicit(g2(x,y), x,-2,2,y,-2,2),
       color=black, key="P",
       point_size=4, point_type = filled_circle,
points([[1/8,-4/5]])
)$

### El método de Lagrange correspondiente a g1

print("Definimos la función de Lagrange auxiliar")$
F1(x,y):=f(x,y)+m*g1(x,y);

Necesitamos calcular los puntos críticos de F1 que pertenecen a g1 

jacobian([F1(x,y)],[x,y]);
addcol(%, [g1(x,y)]);
C1:list_matrix_entries(%);
print("Necesitamos resolver el sistema, para anular las derivadas parciales
y que la solución esté en la gráfica g")$

solve(C1,[x,y,m]);

ev(F1(x,y), x=5/sqrt(509), y=-(sqrt(509)+88)/(4*sqrt(509)), m=(sqrt(509)-40)/40),numer;
ev(F1(x,y), x=-5/sqrt(509), y=-(sqrt(509)-88)/(4*sqrt(509)), m=-(sqrt(509)+40)/40),numer;

wxdraw2d( line_width=3, nticks=10000, proportional_axes=xy,
implicit(g1(x,y), x,-2,2,y,-2,2),
       color=black,
       point_size=4, point_type = filled_circle,
points([[1/8,-4/5]]),
      color=red,
points([[5/sqrt(509),-(sqrt(509)+88)/(4*sqrt(509))]]),
points([[-5/sqrt(509) ,-(sqrt(509)-88)/(4*sqrt(509)) ]])
)$
print("La imagen muestra que corresponden a puntos de mínima y máxima distancia")$

### El método de Lagrange correspondiente a g2

Es análogo