Coordenadas esféricas

kill(all)$
eq1:x-rho*sin(psi)*cos(phi);
eq2:y-rho*sin(psi)*sin(phi);
eq3:z-rho*cos(psi);

jacobian([eq1,eq2,eq3],[rho,psi,phi]);
detsis:determinant(%);

trigsimp(%);

Un ejemplo representativo: F(x,y,u,v)=0

kill(all)$ 
eq1:u+3*x-y^3;
eq2:v+3*y-x^3;

Encontrando u(x,y) y v(x,y)

x0:1$ y0:1$ 
equv1:ev(eq1,x=x0,y=y0);
equv2:ev(eq2,x=x0,y=y0);
sol:solve([equv1,equv2],[u,v]);
u0:ev(u,sol);
v0:ev(v,sol);

[x0,y0,u0,v0];
ev(eq1,x=x0,y=y0,u=u0,v=v0);
ev(eq2,x=x0,y=y0,u=u0,v=v0);

Encontrando la inversa de la función anterior, i.e., x(u,v) e y(u,v)

realonly:true$

u0:-2$ v0:-2$ 
eqxy1:ev(eq1,u=u0,v=v0);
eqxy2:ev(eq2,u=u0,v=v0);
eliminate([eqxy1,eqxy2],[y]);
eliminate([eqxy1,eqxy2],[x]);
solxy:solve([eqxy1,eqxy2],[x,y]);
length(%);

x0:-2$ y0:-2$ 
equv1:ev(eq1,x=x0,y=y0);
equv2:ev(eq2,x=x0,y=y0);
sol:solve([equv1,equv2],[u,v]);
u0:ev(u,sol);
v0:ev(v,sol);

[x0,y0,u0,v0];
ev(eq1,x=x0,y=y0,u=u0,v=v0);
ev(eq2,x=x0,y=y0,u=u0,v=v0);

jacobian([eq1,eq2],[x,y]);
detsis:determinant(%),factor;

/* evaluando el determinante en el punto */
ev(detsis,x=x0,y=y0,sol);

ev(detsis,x=1,y=1,sol);

x0:1$ y0:2$
ev(detsis,x=x0,y=y0,sol);
equv1:ev(eq1,x=x0,y=y0);
equv2:ev(eq2,x=x0,y=y0);
sol:solve([equv1,equv2],[u,v]);
u0:ev(u,sol);
v0:ev(v,sol);

u0:5$ v0:-5$ 
eqxy1:ev(eq1,u=u0,v=v0);
eqxy2:ev(eq2,u=u0,v=v0);
eliminate([eqxy1,eqxy2],[y]);
eliminate([eqxy1,eqxy2],[x]);
solxy:solve([eqxy1,eqxy2],[x,y]);
length(%);

vsxy:solxy[5];
vsxy2:solxy[4];

ev([x,y,u0,v0], vsxy);
ev(eq1,u=u0,v=v0,vsxy);
ev(eq2,u=u0,v=v0,vsxy);

ev([x,y,u0,v0], vsxy2);
ev(eq1,u=u0,v=v0,vsxy2);
ev(eq2,u=u0,v=v0,vsxy2);

Encontrando el polinomio de taylor de orden 2 alrededor de los puntos anteriores

depends([x,y],[u,v]);

eq1u:diff(eq1,u)$
eq1us:subst(['diff(x,u,1)=xu,'diff(x,v,1)=xv,'diff(y,u,1)=yu,'diff(y,v,1)=yv],%);
eq2u:diff(eq2,u)$
eq2us:subst(['diff(x,u,1)=xu,'diff(x,v,1)=xv,'diff(y,u,1)=yu,'diff(y,v,1)=yv],%);
eq1v:diff(eq1,v)$
eq1vs:subst(['diff(x,u,1)=xu,'diff(x,v,1)=xv,'diff(y,u,1)=yu,'diff(y,v,1)=yv],%);
eq2v:diff(eq2,v)$
eq2vs:subst(['diff(x,u,1)=xu,'diff(x,v,1)=xv,'diff(y,u,1)=yu,'diff(y,v,1)=yv],%);

solD1:solve([eq1us,eq2us,eq1vs,eq2vs],[xu,yu,xv,yv]);

factor(%);

eq1u2:diff(eq1,u,2)$
eq1u2s:subst(['diff(x,u,1)=xu,'diff(x,v,1)=xv,'diff(y,u,1)=yu,'diff(y,v,1)=yv, 'diff(y,u,2)=yuu, 'diff(x,u,2)=xuu  ],%);
eq2u2:diff(eq2,u,2)$
eq2u2s:subst(['diff(x,u,1)=xu,'diff(x,v,1)=xv,'diff(y,u,1)=yu,'diff(y,v,1)=yv, 'diff(y,u,2)=yuu, 'diff(x,u,2)=xuu  ],%);
eq1v2:diff(eq1,v,2)$
eq1v2s:subst(['diff(x,u,1)=xu,'diff(x,v,1)=xv,'diff(y,u,1)=yu,'diff(y,v,1)=yv, 'diff(y,v,2)=yvv, 'diff(x,v,2)=xvv  ],%);
eq2v2:diff(eq2,v,2)$
eq2v2s:subst(['diff(x,u,1)=xu,'diff(x,v,1)=xv,'diff(y,u,1)=yu,'diff(y,v,1)=yv, 'diff(y,v,2)=yvv, 'diff(x,v,2)=xvv  ],%);

solD2:solve([eq1u2s,eq2u2s,eq1v2s,eq2v2s],[xuu,yuu,xvv,yvv]);

ev(solD2,solD1);

factor(%);

eq1uv2:diff(eq1,u,1,v,1)$
eq1uv2s:subst(['diff(x,u,1)=xu,'diff(x,v,1)=xv,'diff(y,u,1)=yu,'diff(y,v,1)=yv, 'diff(y,v,1,u,1)=yuv, 'diff(x,v,1,u,1)=xuv  ],%);
eq2uv2:diff(eq2,u,1,v,1)$
eq2uv2s:subst(['diff(x,u,1)=xu,'diff(x,v,1)=xv,'diff(y,u,1)=yu,'diff(y,v,1)=yv, 'diff(y,v,1,u,1)=yuv, 'diff(x,v,1,u,1)=xuv  ],%);

solve([eq2uv2s,eq1uv2s],[xuv,yuv]);
ev(%,solD1),factor;

solD2:solve([eq1u2s,eq2u2s,eq1v2s,eq2v2s,eq2uv2s,eq1uv2s],[xuu,yuu,xvv,yvv,xuv,yuv]);

vsolD2:ev(%,solD1),factor;

Dx: matrix( [xu,xv])$
Dy: matrix( [yu,yv])$

vDx:ev(Dx,solD1);
vDy:ev(Dy,solD1);

D2x: matrix(
 [xuu,xuv], 
 [xuv,xvv]
)$
vD2x:ev(%,vsolD2);

D2y: matrix(
 [yuu,yuv], 
 [yuv,yvv]
)$
vD2y:ev(%,vsolD2);

polx:x0 + vDx.[(u-u0),(v-v0)]+ [(u-u0),(v-v0)].(1/2*vD2x).[(u-u0),(u-v0)];

poly:y0 + vDy.[(u-u0),(v-v0)]+ [(u-u0),(v-v0)].(1/2*vD2y).[(u-u0),(u-v0)];

vsxy;

ev(polx,vsxy);

ev(vDx,vsxy);

ev(vD2x,vsxy);

vsxy2;

ev(polx,vsxy2);

ev(vDx,vsxy2);

ev(vD2x,vsxy2);