# Cálculo diferencial en varias variables

## Mini intro

### Operaciones algebraicas y funciones 

2+2; 3-3; 2*3; 5/10; 3^3;

e;

float(e);

 %e;

float(%e);

 (2+3^2)^3*(5+2^2);

 x:123; y:321; x*y; x/y; x-y;

z=2;

z;

x; y;

x*y;

x+z;

 f(x):= x^2 -x + 1;

f(%pi);

 float(f(%pi));

 ??float;

??renato;

/* MANUAL ONLINE y libro */;

 wxplot2d([f(x)], [x,-5,5]);

wxplot2d([f('x)], ['x,-5,5]);

kill(x);

 wxplot2d([f(x)], [x,-5,5]);

plot2d([f(x),2*sin(x),atan(x)], [x,-2,2])$

log(10);

float(%);

log(%e);

 15!;  
factor(%);

### Resolviendo ecuaciones

kill(all)$

 solve(x^2-a-4,x);

solve(x^2-a-4,a);

 solve(x^2=-1);

 solve(x^5=1,x);

 solve(x^n=1,x);

kk:3$ /* probar con 4,5,6,7, ....*/
solve(x^kk=1,x); 

### Límites de una variable

sin(a*x)/x;

  limit(sin(a*x)/x,x,0);

 limit(exp(1/(1-x)),x,1);

 limit(exp(1/(1-x)),x,1,plus);

 limit(exp(1/(1-x)),x,1,minus);

 limit(log(x),x,0);

### Dibujando en 3D

kill(f,x,y)$
       f(x,y):= sin(x) + cos(y);
       wxplot3d(f(x,y), [x,-5,5], [y,-5,5])$

load(draw)$

wxdraw3d( proportional_axes = 'xyz, xu_grid = 50,  enhanced3d = true,  
    cbrange = [-3,3],   colorbox = true, 
explicit(f(x,y), x,-5,5,y,-5,5) )$

wxdraw3d( proportional_axes = 'xyz, color=blue, xu_grid = 50,  
explicit(f(x,y), x,-5,5,y,-5,5) )$

wxdraw3d( proportional_axes = 'xyz, color=green, xu_grid = 50,  
             explicit(f(x,y), x,-5,5,y,-5,5),
             point_type=filled_circle, point_size =2,color = dark_green, 
             points([[2,1,f(2,1)]])    )$

draw3d( proportional_axes = 'xyz, color=green, xu_grid = 50,  /* view=[10,90], */
             explicit(f(x,y), x,-5,5,y,-5,5),
             point_type=filled_circle, point_size =2,color = dark_green, 
             points([[2,1,f(2,1)]])    )$

## Introducción al cálculo de varias variables

### Continuidad, derivadas parciales y diferenciabilidad

### Ejemplo 1

kill(f)$
f(x,y):=if x^2+y^2=0 then 0 else (x*y)/(x^2+y^2);

limit(f(x,y),x,0); 
limit(f(x,y),y,0);

limit(f(x,m*x),x,0); 

ratsimp((f(x,0)-f(0,0))/x);
ratsimp((f(0,y)-f(0,0))/x);

draw3d( color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]), 
    xu_grid = 100,  enhanced3d = false,  cbrange = [-3,3],    proportional_axes = 'xyz,  
    colorbox = true,  color = blue, explicit(f(x,y), x,-5,5,y,-5,5), 
    color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]))$

draw3d( color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]), 
    xu_grid = 100,  enhanced3d = false,  cbrange = [-3,3],   /* proportional_axes = 'xyz,  */
    colorbox = true,  color = blue, explicit(f(x,y), x,-5,5,y,-5,5), 
    color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]))$

wxdraw2d(color = green, point_size = 5, point_type = filled_circle, points([[0,0]]), 
    yrange=[-.5,1], grid=true, line_width=5,
    color = red,explicit(f(x,x), x,-5,5), 
    color = blue, explicit(f(x,x^2), x,-5,5),
    color=black, explicit(f(x,x^3), x,-5,5))$

### Ejemplo 2

kill(f)$
f(x,y):=if x^2+y^2=0 then 0 else (x^2*y)/(x^2+y^2);

f(1,2); f(0,0);

ev(f(x,y),y=m*x^(1)); /* potencisa x, x^2, x^3 x^(2/3) */
ratsimp(%);radcan(%);

diff(f(x,y), x);

factor(%);

ev(%,y=m*x);
ratsimp(%);

diff(f(x,y), y), factor;

ev(%,y=m*x);
ratsimp(%);

diff(f(x,y));

ratsimp((f(x,0)-f(0,0))/x);
ratsimp((f(0,y)-f(0,0))/x);

/* Diferenciabilidad en (0,0)*/
f(x,y)/sqrt(x^2+y^2);

ev(%,y=m*x),ratsimp;

draw3d( color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]), 
    proportional_axes = 'xyz, xu_grid = 100,  enhanced3d = false,  cbrange = [-3,3],   
    colorbox = true,  color = blue, explicit(f(x,y), x,-5,5,y,-5,5), 
    color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]))$

wxdraw2d(color = green, point_size = 5, point_type = filled_circle, points([[0,0]]), 
    yrange=[-.5,1], grid=true, line_width=5,
    color = red,explicit(f(x,x), x,-5,5), 
    color = blue, explicit(f(x,x^2), x,-5,5),
    color=black, explicit(f(x,x^3), x,-5,5))$

draw3d( color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]), 
    xu_grid = 100,  enhanced3d = false,  cbrange = [-3,3],   proportional_axes = 'xyz, 
    colorbox = true, color = blue,  explicit(f(x,y)/sqrt(x^2+y^2), x,-5,5,y,-5,5),
     color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]))$

draw3d( color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]), 
    xu_grid = 100,  enhanced3d = false,  cbrange = [-3,3],  /* proportional_axes = 'xyz, */
    colorbox = true, color = blue,  explicit(f(x,y)/sqrt(x^2+y^2), x,-5,5,y,-5,5),
     color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]))$

wxdraw2d(color = green, point_size = 5, point_type = filled_circle, points([[0,0]]), 
    yrange=[-.5,1], grid=true, line_width=5,
    color = red,explicit(f(x,x)/sqrt(x^2+x^2), x,-5,5), 
    color = blue, explicit(f(x,x^2)/sqrt(x^2+x^4), x,-5,5),
    color=black, explicit(f(x,x^3)/sqrt(x^2+x^6), x,-5,5))$

f(x,x)/sqrt(x^2+x^2);

### Ejemplo 3

kill(f)$
f(x,y):=if x^2+y^2=0 then 0 else (x^2*y^2)/(x^2+y^2);

f(0,0); f(1,2);

diff(f(x,y), x);

factor(%);

ev(%,y=m*x), ratsimp;

diff(f(x,y), y), factor;

ev(%,y=m*x), ratsimp;

diff(f(x,y));

/* Diferenciabilidad en (0,0)*/
f(x,y)/sqrt(x^2+y^2);

ev(%,y=m*x^2), ratsimp;

radcan(%);

draw3d( color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]), 
    proportional_axes = 'xyz, xu_grid = 100,  enhanced3d = false,  cbrange = [-3,3],   
    colorbox = true,  color = blue, explicit(f(x,y), x,-5,5,y,-5,5), 
    color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]))$

wxdraw2d(color = green, point_size = 5, point_type = filled_circle, points([[0,0]]), 
    yrange=[-.5,1], grid=true, line_width=5,
    color = red,explicit(f(x,x), x,-5,5), 
    color = blue, explicit(f(x,x^2), x,-5,5),
    color=black, explicit(f(x,x^3), x,-5,5))$

draw3d( color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]), 
    proportional_axes = 'xyz, xu_grid = 100,  enhanced3d = false,  cbrange = [-3,3],   
    colorbox = true, color = blue,  explicit(f(x,y)/sqrt(x^2+y^2), x,-5,5,y,-5,5),
     color = black, point_size = 3, point_type = filled_circle, points([[0,0,0]]))$

wxdraw2d(color = green, point_size = 5, point_type = filled_circle, points([[0,0]]), 
    yrange=[-.5,1], grid=true, line_width=5,
    color = red,explicit(f(x,x)/sqrt(x^2+x^2), x,-5,5), 
    color = blue, explicit(f(x,x^2)/sqrt(x^2+x^4), x,-5,5),
    color=black, explicit(f(x,x^3)/sqrt(x^2+x^6), x,-5,5))$

### Ejemplo 4

kill(g)$
g(x,y):=if x^2+y^2=0 then 0 else (x^2+y^2)*sin(1/sqrt(x^2+y^2));

g(1,1);g(0,0);

a:0$ b:0$

g(x,y);

draw3d( enhanced3d = false, xu_grid = 50, color =dark_salmon,
explicit(g(x,y), x,a-1/2,a+1/2,y,b-1/2,b+1/2), 
point_type=filled_circle, point_size =2,color = black, points([[a,b,0]]) )$

draw3d( enhanced3d = false, xu_grid = 50, color = dark_salmon,
explicit(g(x,y), x,a-.2,a+.2,y,b-.2,b+.2), 
point_type=filled_circle, point_size =2,color = black, points([[a,b,0]]) )$

draw3d( enhanced3d = false, xu_grid =100, color = dark_salmon,
explicit(g(x,y), x,a-.1,a+.1,y,b-.1,b+.1), 
point_type=filled_circle, point_size =2,color = black, points([[a,b,0]]) )$

wxplot3d (50*g(x,y), [x,a-.1,a+.1], [y,a-.1,a+.1], [zlabel,""],
[mesh_lines_color,false], [elevation,0], [azimuth,0],
color_bar, [grid,80,80], [ztics,false], [color_bar_tics,1])$

diff(g(x,y),x)$
gx:expand(%);

diff(g(x,y),y)$
gy:expand(%); 
ev(%,y=m*x);

g(x,0);
g(0,y);

/* Derivadas parciales el (0,0) */
gx00:limit((g(x,0)-g(0,0))/x,x,0);
gy00:limit((g(0,y)-g(0,0))/y,y,0);

g(x,y)/sqrt(x^2+y^2);

/* Dibujamos la función (Diferenciabilidad) */
draw3d( enhanced3d = false, xu_grid = 50, color = dark_salmon,
explicit(g(x,y)/sqrt(x^2+y^2), x,a-2,a+2,y,b-2,b+2), 
point_type=filled_circle, point_size = 2,color = blue, points([[a,b,0]]) )$

/* Dibujamos la función (Diferenciabilidad) zoom */
draw3d( enhanced3d = false, xu_grid = 50, color = dark_salmon,
explicit(g(x,y)/sqrt(x^2+y^2), x,a-.2,a+.2,y,b-.2,b+.2), 
point_type=filled_circle, point_size =1,color = blue, points([[a,b,0]]) )$

wxplot3d (50*g(x,y), [x,a-.1,a+.1], [y,a-.1,a+.1], [zlabel,""],
[mesh_lines_color,false], [elevation,0], [azimuth,0],
color_bar, [grid,80,80], [ztics,false], [color_bar_tics,1])$

wxplot3d (10*g(x,y)/sqrt(x^2+y^2), [x,a-.1,a+.1], [y,a-.1,a+.1], [zlabel,""],
[mesh_lines_color,false], [elevation,0], [azimuth,0],
color_bar, [grid,80,80], [ztics,false], [color_bar_tics,1])$

gx;  /* No es continua en (0,0)*/

gy; /* No es continua en (0,0)*/

## Matrices Jacobianas

### Ejemplo 5

kill(f1,f2)$
f1(x,y):= atan( (x+y)/(1-x*y) );
f2(x,y):= tan(x^2/y);

jacobian([f1(x,y),f2(x,y)],[x,y]);

define(jac(x,y),ratsimp(%))$

jac(x,y);

jac(%pi^(1/2),1);

### Ejemplo 6

kill(all)$
f1(x,y,z):= log(1+x^2+y^2+2*z^4);
f2(x,y,z):= tanh(x^2-y*z);

jacobian([f1(x,y,z),f2(x,y,z)],[x,y,z]);

define(jac(x,y,z),ratsimp(%));

jac(1,0,-1);

jac(1,1,1);

define(jac1(x,y,z),jacobian([f1(x,y,z)],[x,y,z]));
define(hes1(x,y,z),hessian(f1(x,y,z),[x,y,z]));
define(jac2(x,y,z),jacobian([f2(x,y,z)],[x,y,z]));
define(hes2(x,y,z),hessian(f2(x,y,z),[x,y,z]));

a:1$ b:1$ c:1$
jac1(a,b,c); hes1(a,b,c); 
jac2(a,b,c); hes2(a,b,c); 

p1:f1(a,b,c)+jac1(a,b,c).[x-a,y-b,z-c]+(1/2)*[x-a,y-b,z-c].hes1(a,b,c).[x-a,y-b,z-c];

p2:f2(a,b,c)+jac2(a,b,c).[x-a,y-b,z-c]+(1/2)*[x-a,y-b,z-c].hes2(a,b,c).[x-a,y-b,z-c];

taylor (f1(x,y,z) , [x, y,z], [a,b,c], 2);
%-p1;

taylor (f2(x,y,z) , [x, y,z], [a,b,c], 2);
%-p2;

## Calculando derivadas direccionales y planos tangentes

Derivadas direccionales

kill(all)$

/* Definimos la norma en R^2 y R^3 */
define(norm2(vv),sqrt(sum(vv[k]^2,k,1,2)));
define(norm3(vv),sqrt(sum(vv[k]^2,k,1,3)));

### Ejemplo 1

/* Definimos la funcion */  
define(f(x,y),x^2*y^3);

/* elegimos el punto */
a:1$ b:2$ c:f(a,b)$

/* Calculamos las derivadas parciales */
diff(f(x,y),x);
diff(f(x,y),y);

/* Calculamos el Jacobiano de la función */
gra:jacobian([f(x,y)],[x,y]);

gran:ev(gra,x=a,y=b);

/* Definimos el vector dirección normalizado */
v1:5$ v2:1$ v:[v1,v2]$  
vec:v/norm2(v);

/* Calculamos la derivada direccional de f segun v y la evaluamos en (a,b) */
Du:gran.vec;  /*   Multiplicación de matrices: probar gra*vec   */

/* Dirección de maxima variación normalizada*/
xgra:gran[1][1]$ ygra:gran[1][2]$
vgra:[xgra,ygra];
ngra:sqrt(  xgra^2+ygra^2 )$
vgra/ngra;

/* Dibujamos la función f(x,y)  y f(x(t),y(t)) en la dirección del vector */
draw3d( enhanced3d = false, border = false, xu_grid = 50, view=[81,0],color = spring_green,  
    explicit(f(x,y), x,a-2,a+2,y,b-2,b+2), /* dibujamos la función */
line_width = 3,  color = dark_green,   /* dibujamos la función según en la dirección del vector */
    parametric(a+v1*t/norm2(v),b+v2*t/norm2(v),f(a+v1*t/norm2(v),b+v2*t/norm2(v)),t,0,2), 
point_type=filled_circle, point_size = 2,color = forest_green, /*punto A*/
    points([[a,b,f(a,b)]]), 
head_angle  = 20, head_length = 0.1, /*  vector  */
    vector([a,b,c],[v1/norm2(v),v2/norm2(v),0]) )$

/* Dibujamos la función f(x,y)  y f(x(t),y(t)) en la dirección del vector y la del gradiente */
draw3d( enhanced3d = false, border = false, xu_grid = 50, view=[81,0], color = spring_green,  
    explicit(f(x,y), x,a-2,a+2,y,b-2,b+2),  /* dibujamos la función */
line_width = 3,  color = dark_green,  /* dibujamos la función según en la dirección del vector */
     parametric(a+v1*t/norm2(v),b+v2*t/norm2(v),f(a+v1*t/norm2(v),b+v2*t/norm2(v)),t,0,2), 
point_type=filled_circle, point_size = 2,color = forest_green, /*punto A*/
    points([[a,b,f(a,b)]]),  head_angle  = 20, head_length = 0.1, /*  vector  */
    vector([a,b,c],[v1/norm2(v),v2/norm2(v),0]),  color = red, /* vector gradiente  */
   vector([a,b,c],[xgra/ngra,ygra/ngra,0]),
     /* dibujamos la función según en la dirección del vector gradiente */
  parametric(a+xgra/ngra*t,b+ygra/ngra*t,f(a+xgra/ngra*t,b+ygra/ngra*t),t,0,2)
)$

### Ejemplo 2

/* EJEMPLO 3: Definimos la funcion */  
kill(f,a,b,c,v1,v2)$
define(f(x,y),sin( (x^2+y^2)/(1+x^2+y^2)) );

/* elegimos el punto */
a:-1$ b:1$ c:f(a,b)$

/* Calculamos las derivadas parciales */
diff(f(x,y),x);
diff(f(x,y),y);

/* Calculamos el Jacobiano de la función */
gra:jacobian([f(x,y)],[x,y]);

gran:ev(gra,x=a,y=b);

/* Definimos el vector dirección normalizado */
v1:1$ v2:3$ v:[v1,v2]$  
vec:v/norm2(v);

/* Calculamos la derivada direccional de f segun v y la evaluamos en (a,b) */
Du:gran.vec;  /*   Multiplicación de matrices: probar gra*vec   */

/* Dirección de maxima variación normalizada*/
xgra:gran[1][1]$ ygra:gran[1][2]$
vgra:[xgra,ygra];
ngra:sqrt(  xgra^2+ygra^2 )$
vgra/ngra;

/* Dibujamos la función */
draw3d( enhanced3d = false, border = false, xu_grid = 50, view=[70,14],color = spring_green,  
    explicit(f(x,y), x,a-2,a+2,y,b-2,b+2), /* dibujamos la función */
line_width = 3,  color = dark_green,   /* dibujamos la función según en la dirección del vector */
    parametric(a+v1*t/norm2(v),b+v2*t/norm2(v),f(a+v1*t/norm2(v),b+v2*t/norm2(v)),t,0,2), 
point_type=filled_circle, point_size = 2,color = forest_green, /*punto A*/
    points([[a,b,f(a,b)]]),
head_angle  = 20, head_length = 0.1, /*  vector  */
    vector([a,b,c],[v1/norm2(v),v2/norm2(v),0]) )$

/* Dibujamos la función f(x,y)  y f(x(t),y(t)) en la dirección del vector y la del gradiente */
draw3d( enhanced3d = false, border = false, xu_grid = 50, view=[70,14], color = spring_green,  
    explicit(f(x,y), x,a-2,a+2,y,b-2,b+2),  /* dibujamos la función */
line_width = 3,  color = dark_green,  /* dibujamos la función según en la dirección del vector */
     parametric(a+v1*t/norm2(v),b+v2*t/norm2(v),f(a+v1*t/norm2(v),b+v2*t/norm2(v)),t,0,2), 
point_type=filled_circle, point_size = 2,color = forest_green, /*punto A*/
    points([[a,b,f(a,b)]]), head_angle  = 20, head_length = 0.1, /*  vector  */
    vector([a,b,c],[v1/norm2(v),v2/norm2(v),0]), color = red, /* vector gradiente  */
   vector([a,b,c],[xgra/ngra,ygra/ngra,0]),
     /* dibujamos la función según en la dirección del vector gradiente */
  parametric(a+xgra/ngra*t,b+ygra/ngra*t,f(a+xgra/ngra*t,b+ygra/ngra*t),t,0,2)
)$

/* Dibujamos la función f(x,y)  y f(x(t),y(t)) en la dirección del vector y la del gradiente y - gradiente*/
draw3d( enhanced3d = false, border = false, xu_grid = 50, view=[70,14], color = spring_green,  
    explicit(f(x,y), x,a-2,a+2,y,b-2,b+2),  /* dibujamos la función */
line_width = 3,  color = dark_green,  /* dibujamos la función según en la dirección del vector */
     parametric(a+v1*t/norm2(v),b+v2*t/norm2(v),f(a+v1*t/norm2(v),b+v2*t/norm2(v)),t,0,2), 
point_type=filled_circle, point_size = 2,color = forest_green, /*punto A*/
    points([[a,b,f(a,b)]]), head_angle  = 20, head_length = 0.1, /*  vector  */
    vector([a,b,c],[v1/norm2(v),v2/norm2(v),0]),color = red, /* vector gradiente  */
   vector([a,b,c],[xgra/ngra,ygra/ngra,0]),
     /* dibujamos la función según en la dirección del vector gradiente */
  parametric(a+xgra/ngra*t,b+ygra/ngra*t,f(a+xgra/ngra*t,b+ygra/ngra*t),t,0,2),
    color = blue, /* vector - gradiente  */   vector([a,b,c],[-xgra/ngra,-ygra/ngra,0]),
     /* dibujamos la función según en la dirección del vector gradiente */
  parametric(a-xgra/ngra*t,b-ygra/ngra*t,f(a-xgra/ngra*t,b-ygra/ngra*t),t,0,1)
)$

### Ejemplo 3

/* Definimos la funcion  */
kill(f,a,b,c)$
define(f(x,y,z),(x/y)^z);

gra:jacobian([f(x,y,z)],[x,y,z]);

v:[2,1,-1]$ 
vec:v/norm3(v);

a:1$ b:1$ c:1$
gran:ev(gra,x=a,y=b,z=c);

Du:gran.vec;

### Ejemplo 4

/*   Definimos la función */
kill(f,a,b,c)$
define(f(x,y,z),(x^2+y^2)^(z));

gra:jacobian([f(x,y,z)],[x,y,z]);

v:[1,-1,1]$ 
vec:v/norm3(v);

a:0$ b:1$ c:1$
gran:ev(gra,x=a,y=b,z=c);

Du:gran.vec;

Planos tangente

kill(all)$
load(draw)$

### Ejemplo 5

kill(all)$
es:1/10$  /* Se reescala por comodidad*/
define(f(x,y),es*(1-(x^2+y^2)));

define(dxf(x,y),diff(f(x,y),x));
define(dyf(x,y),diff(f(x,y),y));

/* Ecuacion del plano */
a:-1$ b:1$ c:f(a,b)$  xmin:a-1$ xmax:a+1$ ymin:b-1$ ymax:b+1$
vn:[dxf(a,b),dyf(a,b),-1]*(-1);
z=f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b);
define(plano(x,y), f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b));

draw3d( proportional_axes = 'xyz, color=green, xu_grid = 100,  /*quitar proportional_axes = 'xyz */
explicit(f(x,y),x,xmin-1,xmax+1,y,ymin-1,ymax+1), color=blue, 
parametric_surface(x,y,plano(x,y),x,xmin,xmax,y,ymin,ymax),
point_type=filled_circle, point_size =2,color = black, points([[a,b,f(a,b)]]), 
head_angle  = 20,head_length = 0.1,vector([a,b,c],vn) )$

### Ejemplo 6

kill(all)$
es:1$
define(f(x,y),es*(  x^2*y^2/(x^2+y^2)  ));

define(dxf(x,y),diff(f(x,y),x));
define(dyf(x,y),diff(f(x,y),y));

/* Ecuacion del plano */
a:-2$ b:1$ c:f(a,b)$  xmin:a-1$ xmax:a+1$ ymin:b-1$ ymax:b+1$
vn:[dxf(a,b),dyf(a,b),-1]*(-1);
z=f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b);
define(plano(x,y), f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b));

draw3d( proportional_axes = 'xyz, color=blue, xu_grid = 25, 
    unit_vectors = true, view=[90,60],
 explicit(f(x,y),x,xmin-1,xmax+1,y,ymin-1,ymax+1), color=cyan, 
parametric_surface(x,y,plano(x,y),x,xmin,xmax,y,ymin,ymax),
point_type=filled_circle, point_size =1,color = black, 
points([[a,b,f(a,b)]]), head_angle  = 20,head_length = 0.1,
vector([a,b,c],vn)  )$

draw3d(color=blue, xu_grid = 25, /*  proportional_axes = 'xyz,  */
    unit_vectors = true, view=[90,60],  title="vector normal?",
 explicit(f(x,y),x,xmin-1,xmax+1,y,ymin-1,ymax+1), color=cyan, 
parametric_surface(x,y,plano(x,y),x,xmin,xmax,y,ymin,ymax),
point_type=filled_circle, point_size =1,color = black, font      = "Arial", 
font_size = 20, points([[a,b,f(a,b)]]), head_angle  = 20,
head_length = 0.1, vector([a,b,c],vn)  )$

### Ejemplo 7

kill(all)$
es:1$
define(f(x,y),es*( atan(x^2+y^2) ));

define(dxf(x,y),diff(f(x,y),x));
define(dyf(x,y),diff(f(x,y),y));

/* Ecuacion del plano */
a:0$ b:1$ c:f(a,b)$  xmin:a-1$ xmax:a+1$ ymin:b-1$ ymax:b+1$
vn:[dxf(a,b),dyf(a,b),-1]*(1);
z=f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b);
define(plano(x,y), f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b));

draw3d( proportional_axes = 'xyz, color=blue, xu_grid = 30,  
      unit_vectors = true, view=[100,110],
explicit(f(x,y),x,xmin-1,xmax+1,y,ymin-1,ymax+1),
color=red, 
parametric_surface(x,y,plano(x,y),x,xmin,xmax,y,ymin,ymax),
point_type=filled_circle, point_size =1,color = black, 
points([[a,b,f(a,b)]]),
head_angle  = 20,head_length = 0.1,vector([a,b,c],vn)
)$

### Ejemplo 8

kill(all)$
es:1$
define(f(x,y),es*((x^2+y^2)*sin((x^2+y^2))  ));

define(dxf(x,y),diff(f(x,y),x));
define(dyf(x,y),diff(f(x,y),y));

/* Ecuacion del plano */
a:0$ b:0$ c:f(a,b)$  xmin:a-1$ xmax:a+1$ ymin:b-1$ ymax:b+1$
vn:[dxf(a,b),dyf(a,b),-1]*(1);
z=f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b);
define(plano(x,y), f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b));

draw3d( proportional_axes = 'xyz, color=blue, xu_grid = 50,  
      unit_vectors = true, 
explicit(f(x,y),x,xmin-.1,xmax+.1,y,ymin-.1,ymax+.1),
color=red, 
parametric_surface(x,y,plano(x,y),x,xmin,xmax,y,ymin,ymax),
point_type=filled_circle, point_size =1,color = black, 
points([[a,b,f(a,b)]]),
head_angle  = 20,head_length = 0.1,vector([a,b,c],-vn) /* +vn */
)$

a:0$ b:1/2$ c:f(a,b)$  xmin:a-1$ xmax:a+1$ ymin:b-1$ ymax:b+1$
vn:[dxf(a,b),dyf(a,b),-1]*(1);
z=f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b);
define(plano(x,y), f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b));

draw3d(color=blue, xu_grid = 50,  proportional_axes = 'xyz,  axis_3d = false,
explicit(f(x,y),x,xmin-.1,xmax+.1,y,ymin-.1,ymax+.1),
unit_vectors = true,  color=red,  zrange = [-1,2.5],  xrange=[-1,1], yrange=[-1,1],
parametric_surface(x,y,plano(x,y),x,xmin,xmax,y,ymin,ymax),
point_type=filled_circle, point_size =1,color = black, points([[a,b,f(a,b)]]),
head_angle  = 20,head_length = 0.1,vector([a,b,c],10*vn)  )$

### Ejemplo 9

kill(all)$
es:2$
define(f(x,y),es*(  exp(-x^2-y^2)  ));

define(dxf(x,y),diff(f(x,y),x));
define(dyf(x,y),diff(f(x,y),y));

/* Ecuacion del plano */
a:0$ b:0$ c:f(a,b)$  xmin:a-1$ xmax:a+1$ ymin:b-1$ ymax:b+1$
vn:[dxf(a,b),dyf(a,b),-1]*(-1);
z=f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b);
define(plano(x,y), f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b));

draw3d( proportional_axes = 'xyz, color=blue, xu_grid = 50,   axis_3d = false,
explicit(f(x,y),x,xmin-1,xmax+1,y,ymin-1,ymax+1), color=cyan, 
parametric_surface(x,y,plano(x,y),x,xmin,xmax,y,ymin,ymax),
point_type=filled_circle, point_size =1,color = black, points([[a,b,f(a,b)]]),
head_angle  = 20,head_length = 0.1,vector([a,b,c],vn) )$

a:1/2$ b:-1/2$ c:f(a,b)$  xmin:a-1$ xmax:a+1$ ymin:b-1$ ymax:b+1$
vn:[dxf(a,b),dyf(a,b),-1]*(-1); 
z=f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b);
define(plano(x,y), f(a,b)+dxf(a,b)*(x-a)+dyf(a,b)*(y-b));

draw3d( color=blue, xu_grid = 50,  proportional_axes = 'xyz, axis_3d = false, /*  proportional_axes = 'xyz */
explicit(f(x,y),x,xmin-1,xmax+1,y,ymin-1,ymax+1), color=red, xrange=[-1.5,1.5],  yrange=[-1.5,1.5],
parametric_surface(x,y,plano(x,y),x,xmin,xmax,y,ymin,ymax), 
point_type=filled_circle, point_size =1,color = black,  points([[a,b,f(a,b)]]),
head_angle  = 20,head_length = 0.1,vector([a,b,c],vn) )$

## Derivadas de orden superior

### Ejemplo 1

kill(all)$
g(x,y):=if x^2+y^2=0 then 0 else atan((x^4+y^4)/(x^2+y^2));

g(1,1);g(0,0);

a:0$ b:0$ /* punto a estudiar */

### Continuidad

g(x,y);

g(x,m*x), ratsimp;

limit(g(x,m*x),x,0);

g(x,m*x^2), ratsimp;

limit(g(x,m*x^2),x,0);

g(x,m*x^(1/2)), ratsimp;
limit(%,x,0);

/*  Estudiando con cuidado el argumento  */
(x^4+y^4)/(x^2+y^2);expand(%);

draw3d( enhanced3d = false, xu_grid = 50, color = blue,
explicit(g(x,y), x,-1,1,y,-1,1), 
point_type=filled_circle, point_size = 3,color = black, points([[a,b,0]]) )$

### Diferenciabilidad

diff(g(x,y),x); 
gx:factor(%);

diff(g(x,y),y)$ 
gy:factor(%);

dd:y^8+2*x^4*y^4+y^4+2*x^2*y^2+x^8+x^4;

aa:y^8+2*x^4*y^4+x^8,factor;

bb:y^4+2*x^2*y^2+x^4,factor;

aa+bb-dd,ratsimp;

g(x,0); g(0,y);

/* Derivadas parciales el (0,0) */
gx00:limit((g(x,0)-g(0,0))/x,x,0);
gy00:limit((g(0,y)-g(0,0))/y,y,0);

/* Dibujamos la función (Diferenciabilidad) */
draw3d( enhanced3d = false, xu_grid = 50, color = red,
explicit(g(x,y)/sqrt(x^2+y^2), x,-1,1,y,-1,1), 
point_type=filled_circle, point_size = 1,color = black, points([[a,b,0]]) )$

 /* Dibujamos la función   (Diferenciabilidad) zoom */
	draw3d (enhanced3d = false, xu_grid = 50,
	color = red, explicit(g(x,y)/sqrt(x^2+y^2), x,-.1,.1,y,-.1,.1),
	point_type=filled_circle, point_size = 1,color = black,
	points([[0,0,0]]) )$

gd(x,y):=g(x,y)/sqrt(x^2+y^2);

gd(x,m*x),ratsimp;

limit(gd(x,m*x),x,0);

gd(x,m*x^2),ratsimp;

limit(gd(x,m*x^2),x,0);

gd(x,m*sqrt(x) ),ratsimp;

limit( gd(x, m*sqrt(x) ) , x, 0);

gd(x,y);

((y^4+x^4)/(y^2+x^2))/sqrt(y^2+x^2);

expand((y^4+x^4)/(y^2+x^2)^(3/2));

y^4/(y^2+x^2)^(3/2)+x^4/(y^2+x^2)^(3/2);

y^4/(y^2)^(3/2)+x^4/(x^2)^(3/2);

### Derivadas de orden 2

diff(g(x,y),x,2);

dgxx:factor(%);

diff(g(x,y),y,2);
dgyy:factor(%);

diff(g(x,y),x,1,y,1);

dgxy:factor(%);

/* Continuidad de las segundas derivadas */
ev(dgxy,y=m*x), ratsimp;

limit(%,x,0);

ev(dgxx,y=m*x), ratsimp;
limit(%,x,0);

/* Diferenciabilidad de las derivadas parciales */
ev(gx,y=0);
ev(gx,x=0);

(ev(gx,y=0)-gx00)/x;

limit(%,x,0);

(gx-2*x)/sqrt(x^2+y^2);
ratsimp(%)  /* definicion de diferenciabilidad en (0,0)*/;

ev(%,y=m*x);

ratsimp(%);

limit(%,x,0);

ev(gy,y=0);
ev(gy,x=0);

(gy-2*y)/sqrt(x^2+y^2);
ratsimp(%);

ev(%,y=m*x);ratsimp(%);

limit(%,x,0);

### Ejemplo 2

kill(all)$
g(x,y):=if x+y=0 then 0 else x*y*sin((1/2)*(x-y)/(x+y));

g(-1,1);g(0,0);g(2,1);

a:0$ b:0$ /* punto a estudiar ¿Qué pasa en (a-a), a distinto de 0? */

### Continuidad en (0,0). 

g(x,y);

draw3d( enhanced3d = false, xu_grid = 50, color = blue,
explicit(g(x,y), x,-1,1,y,-1,1), 
point_type=filled_circle, point_size = 3,color = black, points([[a,b,0]]) )$

### Diferenciabilidad

gx:diff(g(x,y),x); 

x*y*(1/(2*(y+x))-(x-y)/(2*(y+x)^2));

factor(%);

gy:diff(g(x,y),y);

x*y*(-1/(2*(y+x))-(x-y)/(2*(y+x)^2));
aa:factor(%);

ev(aa,x=1/n,y=-1/n+1/n^2);
ratsimp(%);

g(x,0); g(0,y);

/* Derivadas parciales el (0,0) */
gx00:limit((g(x,0)-g(0,0))/x,x,0);
gy00:limit((g(0,y)-g(0,0))/y,y,0);

/* Dibujamos la función (Diferenciabilidad) */
draw3d( enhanced3d = false, xu_grid = 50, color = red,
explicit(g(x,y)/sqrt(x^2+y^2), x,-1,1,y,-1,1), 
point_type=filled_circle, point_size = 1,color = black, points([[a,b,0]]) )$

 /* Dibujamos la función   (Diferenciabilidad) zoom */
	draw3d (enhanced3d = false, xu_grid = 50,
	color = red, explicit(g(x,y)/sqrt(x^2+y^2), x,-.1,.1,y,-.1,.1),
	point_type=filled_circle, point_size = 1,color = black,
	points([[0,0,0]]) )$

gd(x,y):=g(x,y)/sqrt(x^2+y^2);

gd(x,y),ratsimp;

(x*y)/sqrt(y^2+x^2); /* ¿tiene límite si (x,y)→(0,0)? */

### Derivadas de orden 2

diff(g(x,y),x,2);

dgxx:ratsimp(%);

diff(g(x,y),y,2)$
dgyy:factor(%);

diff(g(x,y),x,1,y,1);

dgxy:factor(%);

/* Continuidad de las segundas derivadas */
ev(dgxy,y=m*x), ratsimp;

limit(%,x,0), ratsimp;

ev(dgxx,y=m*x)$;
limit(%,x,0), ratsimp;

/* Calculo de las derivadas parciales en (0,0) */
ev(gx,y=0);
ev(gx,x=0);

ev(gy,y=0);
ev(gy,x=0);

(ev(gy,y=0)-gy00)/x;
gxy00:limit(%,x,0);

(ev(gx,x=0)-gx00)/y;
gyx00:limit(%,y,0);

(ev(gx,y=0)-gx00)/x;
gxx00:limit(%,x,0);

(ev(gy,x=0)-gy00)/y;
gyy00:limit(%,x,0);

### Ejemplo 3

kill(all)$
g(x,y):=if x*y=0 then 0 else (sin(x*y)-x*y)/(x^2*y);

g(0,1);g(0,0);g(1,0);g(1,1);

a:0$ b:0$ /* punto a estudiar ¿Qué pasa en (a,0) y (0,a), a distinto de 0? */

### Continuidad en (0,0). 

g(x,y);

draw3d( enhanced3d = false, xu_grid = 50, color = blue,
explicit(g(x,y), x,-1,1,y,-1,1), 
point_type=filled_circle, point_size = 3,color = black, points([[a,b,0]]) )$

taylor(sin(x*y),[x,y],[0,0],7);

taylor(g(x,y),[x,y],[0,0],7);

-(y^3*x^3)/6/(x^2*y);

### Diferenciabilidad

gx:diff(g(x,y),x); 

ev(gx, [sin(x*y)=taylor(sin(x*y),[x,y],[0,0],7), cos(x*y)=taylor(cos(x*y),[x,y],[0,0],7)]);

gy:diff(g(x,y),y);

ev(gy, [sin(x*y)=taylor(sin(x*y),[x,y],[0,0],7), cos(x*y)=taylor(cos(x*y),[x,y],[0,0],7)]);

g(x,0); g(0,y);

/* Derivadas parciales el (x_0,0) */
gxx_00:limit((g(x_0+h,0)-g(0,0))/h,h,0);
gyx_00:limit((g(x0,h)-g(0,0))/h,h,0);

/* Derivadas parciales el (0,y_0) */
gx0y_0:limit((g(h,y0)-g(0,0))/h,h,0);
gy0y_0:limit((g(0,y0+h)-g(0,0))/h,h,0);

/* Dibujamos la función (Diferenciabilidad) en (0,0) */
draw3d( enhanced3d = false, xu_grid = 50, color = red,
explicit(g(x,y)/sqrt(x^2+y^2), x,-1,1,y,-1,1), 
point_type=filled_circle, point_size = 1,color = black, points([[a,b,0]]) )$

 /* Dibujamos la función   (Diferenciabilidad eb (0,0)) zoom */
	draw3d (enhanced3d = false, xu_grid = 50,
	color = red, explicit(g(x,y)/sqrt(x^2+y^2), x,-.1,.1,y,-.1,.1),
	point_type=filled_circle, point_size = 1,color = black,
	points([[0,0,0]]) )$

gd(x,y):=g(x,y)/sqrt(x^2+y^2);

gd(x,y),ratsimp;

ev(%, [sin(x*y)=taylor(sin(x*y),[x,y],[0,0],7), cos(x*y)=taylor(cos(x*y),[x,y],[0,0],7)]);

factor(%);

radcan(%);

(x*y^2)/sqrt(y^2+x^2); /* ¿tiene límite si (x,y)→(0,0)? */

### Derivadas de orden 2 en (0,0)

/* Calculo de las derivadas parciales en (0,0) */
ev(gx,y=0);

ev(gx,x=0);

/* Hay que redefinir las primeras derivadas con los valores en (0,a) y (a,0)  con a distinto de 0 */

(ev(gx,y=0)-gx00)/x;

gxx:diff(g(x,y),x,2);

ev(gxx, [sin(x*y)=taylor(sin(x*y),[x,y],[0,0],7), cos(x*y)=taylor(cos(x*y),[x,y],[0,0],7)]);

ev(gxx, [sin(x*y)=taylor(sin(x*y),[x,y],[0,0],9), cos(x*y)=taylor(cos(x*y),[x,y],[0,0],9)]);

gyy:diff(g(x,y),y,2);

ev(gyy, [sin(x*y)=taylor(sin(x*y),[x,y],[0,0],9), cos(x*y)=taylor(cos(x*y),[x,y],[0,0],9)]);

gxy:diff(g(x,y),x,1,y,1);

ev(gxy, [sin(x*y)=taylor(sin(x*y),[x,y],[0,0],9), cos(x*y)=taylor(cos(x*y),[x,y],[0,0],9)]);

taylor(g(x,y),[x,y],[0,0],9);