/*
 * This file:
 *   http://anggtwu.net/MAXIMA/2024-2-C3-P1-Q2.mac.html
 *   http://anggtwu.net/MAXIMA/2024-2-C3-P1-Q2.mac
 *          (find-angg "MAXIMA/2024-2-C3-P1-Q2.mac")
 * Author: Eduardo Ochs <eduardoochs@gmail.com>
 *
 * (defun e () (interactive) (find-angg "MAXIMA/2024-2-C3-P1-Q2.mac"))
*/

* (setq last-kbd-macro (kbd "C-a %M SPC C-a <down>"))
* (eepitch-maxima)
* (eepitch-kill)
* (eepitch-maxima)
load("/usr/share/emacs/site-lisp/maxima/emaxima.lisp")$
load("~/MAXIMA/barematrix1.mac")$
:lisp (setf (get '$display2d 'assign) nil)
display2d:'emaxima$

* (eepitch-maxima)
* (eepitch-kill)
* (eepitch-maxima)
hessian(x^3*y^4, [x,y]);
linel : 110;
linenum : 0;
** (eeks)
* (eeks 0.02 0.02)
gradef(W  (x,y), W_x (x,y), W_y (x,y))$
gradef(W_x(x,y), W_xx(x,y), W_xy(x,y))$
gradef(W_y(x,y), W_xy(x,y), W_yy(x,y))$

dd(F)       := [F,
                diff(F,x),   diff(F,y),
                diff(F,x,2), diff(F,x,1,y,1), diff(F,y,2)]$
aa(o,x0y0)  := at(o, [x=x0y0[1], y=x0y0[2]])$
mm(abcdef)  := block([a,b,c,d,e,f],
                     [a,b,c,d,e,f]:abcdef,
                     [a,b*Dx,c*Dy,d*Dx^2/2,e*Dx*Dy,f*Dy^2/2])$
ss(abcdef)  := block([a,b,c,d,e,f],
                     [a,b,c,d,e,f]:abcdef,
                     a+b+c+d+e+f)$
toM(abcdef) := block([a,b,c,d,e,f],
                     [a,b,c,d,e,f]:abcdef,
                     matrix([a,"",""], [b,c,""], [d,e,f]))$

D2       (F) := toM(dd(F));
D2at(x0y0,F) := toM(aa(dd(F),x0y0));
T2M (x0y0,F) := toM(mm(aa(dd(F),x0y0)));
T2  (x0y0,F) :=  ss(mm(aa(dd(F),x0y0)));

/* Alguns testes: */
          dd(W(x,y));
      toM(dd(W(x,y)));
      toM([1,2,3,4,5,6]);
       aa(dd(W(x,y)),[3,4]);
   toM(aa(dd(W(x,y)),[3,4]));
    mm(aa(dd(W(x,y)),[3,4]));
toM(mm(aa(dd(W(x,y)),[3,4])));
 ss(mm(aa(dd(W(x,y)),[3,4])));

D2        (W(x,y));
D2at([3,4],W(x,y));
T2M ([3,4],W(x,y));
T2  ([3,4],W(x,y));

F  : x*y*(6 - 2*x - y);
F  : expand(F);
P1 : [0,6];
P2 : [1,2];
P3 : [3,0];
P4 : [0,0];

/* (2a) */
D2(W(x,y));
F;
D2F : D2(F);

/* (2b) */
[P1, D2F, D2FP1:D2at(P1,F)];
[P2, D2F, D2FP2:D2at(P2,F)];
[P3, D2F, D2FP3:D2at(P3,F)];
[P4, D2F, D2FP4:D2at(P4,F)];

/* (2c) */
TM2([x0,y0],W(x,y));
[P1, D2FP1, T2M(P1,F), T2(P1,F)];
[P2, D2FP2, T2M(P2,F), T2(P2,F)];
[P3, D2FP3, T2M(P3,F), T2(P3,F)];
[P4, D2FP4, T2M(P4,F), T2(P4,F)];

/* (2d) */
grad(F) := [diff(F,x),diff(F,y)]$
   H(F) := hessian(F, [x,y])$
detH(F) := determinant(H(F))$
crit(F) := [F, grad(F), H(F), detH(F)]$
crit(F) := matrix([F, grad(F)], [H(F), detH(F)])$
   crit(W(x,y));
aa(crit(F), P1);
aa(crit(F), P2);
aa(crit(F), P3);
aa(crit(F), P4);

/* definicao da funcao G */
P2;
T2(P2,F);
G__ :                        T2(P2,F);
G_  : subst([Dx=x-1,Dy=y-2], T2(P2,F));
G   : expand(G_);

/* definicao dos conjuntos B e C */
inC_(x,y)   := y <= 6 - 2*x;
inC (x,y,o) := if inC_(x,y) then o else "";
numsB(expr) :=
  apply(matrix,
        makelist(makelist(ev(expr), x,0,3),
                 y, seqby(6,0,-1)))$
numsC(expr) :=
  apply(matrix,
        makelist(makelist(inC(x,y,ev(expr)), x,0,3),
                 y, seqby(6,0,-1)))$
numsB([x,y]);
numsC([x,y]);

/* (2e) */
F;
[numsC(   x  *y),
 numsC(   x^2*y),
 numsC(   x  *y^2)];
[numsC( 6*x  *y),
 numsC(-2*x^2*y),
 numsC(  -x  *y^2),
 numsC(F)];

/* (2f) */
G__;
G_;
Dx : x-1;
Dy : y-2;
[numsC(Dx^2),
 numsC(Dx*Dy),
 numsC(Dy^2)];
[numsC(4),
 numsC(-4*Dx^2),
 numsC(-2*Dx*Dy),
 numsC(  -Dy^2)];
numsC(G);

/* (2g) */
grad(F);

Fx : diff(F,x);
 -4*x*y + 6*y - y^2;
[numsC(-4*x*y),
 numsC(6*y),
 numsC(-y^2),
 numsC(Fx)];

Fy : diff(F,y);
 -2*x*y + 6*x - 2*x^2;
[numsC(-2*x*y),
 numsC(6*x),
 numsC(-2*x^2),
 numsC(Fy)];

numsC(grad(F));

*



** (c3m242p1p 5 "questao-2" "(2g)")
** (c3m242p1a   "questao-2"   "g)")
** (find-myqdraw "myqdraw3.mac")
** (find-myqdraw "myqdraw3.mac" "myimp")
load_myqdraw()$
[topdf_a,topdf_b,topdf_c,topdf_n] : ["~/LATEX/", "2024-2-C3/", "P1-Q2-v2", 0];
myqdrawp_to_new_pdf(); myps(s):=ps(s/5);
myqdrawp_to_screen (); myps(s):=ps(s);

[xmin,ymin, xmax,ymax] : [-1,-1, 4,7]$
level(zz,color) := myimp1(F=zz, lc(color), lk(z=zz))$
levels() := [level( 3.95, gray),
             level( 3.90, gray),
             level( 3.85, gray),
             level( 3,    red),
             level( 2,    orange),
             level( 1,    gold),
             level( 0,    forest_green),
             level(-1,    blue)]$

** As curvas de nível da F, sem os vetores gradientes:
level(zz,color) := myimp1(F=zz, lc(color), lk(z=zz))$
F_levels : myqdrawp(xyrange(), levels());

** As curvas de nível da G, sem os vetores gradientes:
level(zz,color) := myimp1(G=zz, lc(color), lk(z=zz))$
G_levels : myqdrawp(xyrange(), levels());

** Os conjuntos B e C:
B : create_list([x,y], y,seqby(6,0,-1), x,seq(0,3));
eq1 : y = 6 - 2*x;
eq2 : solve(eq1,x);
subst(eq2, x);
define(xmaxC(y), subst(eq2, x));
C : create_list([x,y], y,seqby(6,0,-1), x,seq(0,xmaxC(y)));

atxy(oo,xy) := aa(oo,xy);

** O gradiente da F nos pontos de C:
** (find-myqdraw "myqdraw3.mac" "myPv")
define(v_at (xy), atxy(grad(F)/10,xy))$
define(Pv_at(xy), myPv(xy,v_at(xy),[myps(1)],hl(0.15)))$
Pv_at([1,3]);
grads_at_C : makelist(Pv_at(xy), xy, C)$
F_grads : myqdrawp(xyrange(), grads_at_C);

** O gradiente da F nos pontos de C e as curvas de nível:
level(zz,color) := myimp1(F=zz, lc(color), lk(z=zz))$
F_all : myqdrawp(xyrange(), levels(), grads_at_C);

matrix([F_levels, G_levels, F_grads, F_all]);