## Ejemplo TFI sistema

eq1:(x-1)^2+y^2-z;
eq2:x^2+y^2+z^2-1;

Haciendo "trampas": resolviendo respecto a x:  ( y(x), z(x) )

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

s1:rhs(solx[1][2]); s2:rhs(solx[3][2]);

num(s1^2);

wxplot2d([ num(s1^2) , num(s2^2) ] ,[x,-1,9/8], [y,-14, 3]);

### Comprobando que hay soluciones del sistema 

ev(solx,x=1)$
radcan(%);

ev(solx,x=1/2)$
help0:radcan(%);

ev(solx,x=5/8)$
help1:radcan(%);

### Calculando las matrices jacobianas y derivadas y comprobando Fy' es invertible

Fy:jacobian([eq1,eq2],[y,z]);
Fx:jacobian([eq1,eq2],[x]);

determinant(Fy);

factor(%);

### Encontrando la expresión de las derivadas 

invert(Fy);

invert(Fy).Fx;

df:ratsimp(%);

help0[1];

ev(df,x=1/2, help0[1] ),ratsimp;

help1[1];

ev(df,x=1/2, help0[1] ),ratsimp;

ev(df,x=5/8, help1[1] ),ratsimp;

### El otro camino

/*A partir de aquí y=y(x) y z=z(x) */
depends(y,x);
depends(z,x);

diff(eq1,x);
deq1:subst(['diff(z,x,1)=zx,'diff(y,x,1)=yx],%);
diff(eq2,x);
deq2:subst(['diff(z,x,1)=zx,'diff(y,x,1)=yx],%);

derf:solve([deq1,deq2],[yx,zx]);

diff(eq1,x,2);
d2eq1:subst(['diff(z,x,1)=zx,'diff(y,x,1)=yx, 'diff(z,x,2)=zxx, 'diff(y,x,2)=yxx],  %);

diff(eq2,x,2);
d2eq2:subst(['diff(z,x,1)=zx,'diff(y,x,1)=yx, 'diff(z,x,2)=zxx, 'diff(y,x,2)=yxx],  %);

der2f:solve([d2eq1,d2eq2],[yxx,zxx]);

derf[1];

subst(derf[1],der2f);

ratsimp(%);

Ejercicio: Una opción mejor: x(z) e y(z)

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

## Ejemplo: y²z+x log(z)-x=0

### Definimos la función y compromabos las hipótesis del TFI

kill(all)$
diff(z,x);diff(z,y);

F(x,y,z):=y^2*z+x*log(z)-x;
a:1$ b:-1$c:1$

F(a,b,c);

diff(F(x,y,z),z);
ev(%,x=a,y=b,z=c);

### Calculamos las primeras derivadas

depends(z,[x,y]);

#### Calculamos zx

eq1:diff(F(x,y,z),x);

solve(eq1,'diff(z,x,1));
zx:rhs(%[1]);
vzx:ev(%,x=a,y=b,z=c);

#### Calculamos zy

eq2:diff(F(x,y,z),y);

solve(eq2,'diff(z,y,1));
zy:rhs(%[1]);
vzy:ev(%,x=a,y=b,z=c);

### Calculamos las segundas derivadas

#### Calculando zxx

zxx1:diff(F(x,y,z),x,2);

eq1*z,ratsimp;  /* Sin denominadores como en la pizarra */
d2eq1:diff(%,x),expand;

subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx, 'diff(z,x,2)=zxx],   d2eq1);

ev(%, x=a,y=b,z=c);

solve(%,zxx);

subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx, 'diff(z,x,2)=zxx],   d2eq1);

solve(%,zxx);

ev(%, x=a,y=b,z=c);

/*  Derivando dos veces */
zxx1:diff(F(x,y,z),x,2);
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx, 'diff(z,x,2)=zxx],   zxx1);
ev(%, x=a,y=b,z=c);
solve(%,zxx);

/*Método general */
eqzxx:solve(zxx1,'diff(z,x,2));

/* calculando zxx  el general */
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx],eqzxx);
rhs(%[1])$
vzxx: ev(%, x=a,y=b,z=c);

/* Usando la fórmula del TFI */
diff(zx,x);
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx],%);

vzxx: ev(%, x=a,y=b,z=c);

/* calculando zxx directamente en el punto (a,b,c) */
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx,x=a,y=b,z=c],rhs(eqzxx[1]));
vzxx:subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx,x=a,y=b,z=c],rhs(eqzxx[1]));

#### Calculando zyy

zyy1:diff(F(x,y,z),y,2);

subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx, 'diff(z,y,2)=zyy],   zyy1);
ev(%, x=a,y=b,z=c);
solve(%,zyy);

eqzyy:solve(diff(F(x,y,z),y,2),'diff(z,y,2));

/* calculando zyy  el general */
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx],eqzyy);
rhs(%[1])$
vzyy: ev(%, x=a,y=b,z=c);

/* calculando zyy en el punto (a,b,c) */
subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx],eqzyy);
vzyy:subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx,x=a,y=b,z=c],rhs(eqzyy[1]));

#### Calculando zxy

zxy1:diff(F(x,y,z),x,1,y,1);

subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx,'diff(z,x,1,y,1)=zxy],   zxy1);
ev(%, x=a,y=b,z=c);
solve(%,zxy);

eqzxy:solve(diff(F(x,y,z),x,1,y,1),'diff(z,x,1,y,1));

subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx],eqzxy);
rhs(%[1])$
vzxy: ev(%, x=a,y=b,z=c);

/* calculando zyy en el punto (a,b,c) */
subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx],eqzxy);
vzxy:subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx,x=a,y=b,z=c],rhs(eqzxy[1]));

### Calculando el polinomio de Taylor

#### La diferencial Dz(x,y)

Dz: matrix( [zx,zy] );

vDz:ev(Dz, x=a,y=b,z=c);

#### La diferencial D²z(x,y)

D2z: matrix(
 [zxx,zxy], 
 [zxy,zyy]
);

vD2z:ev(D2z,x=a,y=b,z=c,zxx=vzxx,zyy=vzyy,zxy=vzxy);

#### El polinomio de orden 2: z(a,b)+ Dz(a,b)(x,y)+(1/2)D²z(a,b)(x,y)

c + vDz.[(x-a),(y-b)]+ [(x-a),(y-b)].(1/2*vD2z).[(x-a),(y-b)];

### Cáclulo derivando las expresiones del TFI

kill(all)$
F(x,y,z):=y^2*z+x*log(z)-x;
a:1$ b:-1$c:1$

define(zx(x,y),-diff(F(x,y,z),x)/diff(F(x,y,z),z));
define(zy(x,y),-diff(F(x,y,z),y)/diff(F(x,y,z),z));

vzx:ev(zx(x,y),x=a,y=b,z=c );
vzy:ev(zy(x,y),x=a,y=b,z=c );

depends(z,[x,y]);

dzxx:diff(zx(x,y),x);
subst(['diff(z,x,1)=zx],dzxx);
ev(%, x=a,y=b,z=c,zx=vzx);

dzyy:diff(zy(x,y),y)$
subst(['diff(z,y,1)=zy],dzyy);
ev(%, x=a,y=b,z=c,zy=vzy);

dzxy:diff(zy(x,y),x)$
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx],dzxy);
ev(%, x=a,y=b,z=c,zx=vzx,zy=vzy);

## Ejemplo: 3 x y z- z³ -1=0

### Definimos la función y compromabos las hipótesis del TFI

kill(all)$
diff(z,x);diff(z,y);

F(x,y,z):=3*x*y*z-z^3+1;
a:0$ b:1$c:1$

F(a,b,c);

diff(F(x,y,z),z);
ev(%,x=a,y=b,z=c);

### Calculamos las primeras derivadas

depends(z,[x,y]);

#### Calculamos zx

eq1:diff(F(x,y,z),x);

solve(eq1,'diff(z,x,1));
zx:rhs(%[1]);
vzx:ev(%,x=a,y=b,z=c);

#### Calculamos zy

eq2:diff(F(x,y,z),y);

solve(eq2,'diff(z,y,1));
zy:rhs(%[1]);
vzy:ev(%,x=a,y=b,z=c);

### Calculamos las segundas derivadas

#### Calculando zxx

zxx1:diff(F(x,y,z),x,2);

eq1*z,ratsimp;  /* Sin denominadores como en la pizarra */
d2eq1:diff(%,x),expand;

subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx, 'diff(z,x,2)=zxx],   d2eq1);

ev(%, x=a,y=b,z=c);

solve(%,zxx);

subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx, 'diff(z,x,2)=zxx],   d2eq1);

solve(%,zxx);

ev(%, x=a,y=b,z=c);

/*  Derivando dos veces */
zxx1:diff(F(x,y,z),x,2);
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx, 'diff(z,x,2)=zxx],   zxx1);
ev(%, x=a,y=b,z=c);
solve(%,zxx);

/*Método general */
eqzxx:solve(zxx1,'diff(z,x,2));

/* calculando zxx  el general */
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx],eqzxx);
rhs(%[1])$
vzxx: ev(%, x=a,y=b,z=c);

/* Usando la fórmula del TFI */
diff(zx,x);
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx],%);

vzxx: ev(%, x=a,y=b,z=c);

/* calculando zxx directamente en el punto (a,b,c) */
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx,x=a,y=b,z=c],rhs(eqzxx[1]));
vzxx:subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx,x=a,y=b,z=c],rhs(eqzxx[1]));

#### Calculando zyy

zyy1:diff(F(x,y,z),y,2);

subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx, 'diff(z,y,2)=zyy],   zyy1);
ev(%, x=a,y=b,z=c);
solve(%,zyy);

eqzyy:solve(diff(F(x,y,z),y,2),'diff(z,y,2));

/* calculando zyy  el general */
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx],eqzyy);
rhs(%[1])$
vzyy: ev(%, x=a,y=b,z=c);

/* calculando zyy en el punto (a,b,c) */
subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx],eqzyy);
vzyy:subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx,x=a,y=b,z=c],rhs(eqzyy[1]));

#### Calculando zxy

zxy1:diff(F(x,y,z),x,1,y,1);

subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx,'diff(z,x,1,y,1)=zxy],   zxy1);
ev(%, x=a,y=b,z=c);
solve(%,zxy);

eqzxy:solve(diff(F(x,y,z),x,1,y,1),'diff(z,x,1,y,1));

subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx],eqzxy);
rhs(%[1])$
vzxy: ev(%, x=a,y=b,z=c);

/* calculando zxy en el punto (a,b,c) */
subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx],eqzxy);
vzxy:subst(['diff(z,y,1)=vzy,'diff(z,x,1)=vzx,x=a,y=b,z=c],rhs(eqzxy[1]));

### Calculando el polinomio de Taylor

#### La diferencial Dz(x,y)

Dz: matrix( [zx,zy] );

vDz:ev(Dz, x=a,y=b,z=c);

#### La diferencial D²z(x,y)

D2z: matrix(
 [zxx,zxy], 
 [zxy,zyy]
);

vD2z:ev(D2z,x=a,y=b,z=c,zxx=vzxx,zyy=vzyy,zxy=vzxy);

#### El polinomio de orden 2: z(a,b)+ Dz(a,b)(x,y)+(1/2)D²z(a,b)(x,y)

c + vDz.[(x-a),(y-b)]+ [(x-a),(y-b)].(1/2*vD2z).[(x-a),(y-b)];

### Cálculo derivando las expresiones del TFI

kill(all)$
F(x,y,z):=3*x*y*z-z^3+1;
a:0$ b:1$c:1$

define(zx(x,y),-diff(F(x,y,z),x)/diff(F(x,y,z),z));
define(zy(x,y),-diff(F(x,y,z),y)/diff(F(x,y,z),z));

vzx:ev(zx(x,y),x=a,y=b,z=c );
vzy:ev(zy(x,y),x=a,y=b,z=c );

depends(z,[x,y]);

dzxx:diff(zx(x,y),x);
subst(['diff(z,x,1)=zx],dzxx);
ev(%, x=a,y=b,z=c,zx=vzx);

dzyy:diff(zy(x,y),y)$
subst(['diff(z,y,1)=zy],dzyy);
ev(%, x=a,y=b,z=c,zy=vzy);

dzxy:diff(zy(x,y),x)$
subst(['diff(z,y,1)=zy,'diff(z,x,1)=zx],dzxy);
ev(%, x=a,y=b,z=c,zx=vzx,zy=vzy);