How to expand a composite function into series?

I need to expand such a function $$g[y,z(x,y)]=\frac{-y (z+1)^4-z^4-4 z^3+8 z+8}{z+1},\tag{1}$$ into powers of $$x$$ and $$y$$. Among $$x,y,z$$ there is a constraint equation, for example $$(3 y+3) z^4+z^3 (4-8 x+12 y)+z^2 (18 y-24 x)+z (12 y-24 x)-8 x+3 y=0.\tag{2}$$ There is also the requirement, $$z\rightarrow 0$$ as $$x,y \rightarrow 0$$.

The most direct method is to solve Eq.(2) for $$z$$, and then substituting the $$z=z(x,y)$$ into Eq.(1) to obtain $$f(x,y)$$. At last, one can use, such as

Series[f[x,y],{x,0,2},{y,0,1}]

to expand the composite function as a series of $$x$$ and $$y$$.

The code for the two equations are:

g[y_,z_]:=(8 + 8 z - 4 z^3 - z^4 - y (1 + z)^4)/(1 + z);

and

-8 x + 3 y + (-24 x + 12 y) z + (-24 x + 18 y) z^2 + (4 - 8 x +
12 y) z^3 + (3 + 3 y) z^4 == 0;

I tried this direct method. It takes a very long time, but I do not get a result. If the constraint equation is not only a quartic equation, but a more complicated algebraic equation, the direct method may not work. I want to know whether there are simpler methods to deal with this kind of problems?

According to @Daniel Lichtblau's method, one can first expand the implicit function equation (the constraint equation) into a series of $$x,y$$. The key problem here is that we cannot first expand $$z=z(x,y)$$ formally as a standard Taylor series, and then substitute it into the constraint equation to solve the coefficients, because the true series of $$z=z(x,y)$$ may have minus power, even fractional power of $$x$$ or $$y$$. Below is an example: a simple quadratic equation,

-x + y + (4 x + y) z - (3 + 2 y) z^2==0

We certainly can solve it to obtain $$z=\frac{\sqrt{16 x^2-12 x+9 y^2+12 y}-4 x-y}{2 (-2 y-3)}.$$ (* the first root *) Expanding the solution, Series[z, {x, 0, 2}, {y, 0, 1}], one can obtain such a result $$\left(-\frac{\sqrt{y}}{\sqrt{3}}+\frac{y}{6}+O\left(y^{3/2}\right)\right)+x \left(\frac{1}{2 \sqrt{3} \sqrt{y}}+\frac{2}{3}-\frac{25 \sqrt{y}}{48 \sqrt{3}}-\frac{4 y}{9}+O\left(y^{3/2}\right)\right)+x^2 \left(\frac{1}{8 \sqrt{3} y^{3/2}}-\frac{19 \sqrt{3}}{64 \sqrt{y}}+\frac{333 \sqrt{3} \sqrt{y}}{1024}+O\left(y^{3/2}\right)\right)+O\left(x^3\right).$$

I also tried the methods on this post: [Calculating Taylor polynomial of an implicit function given by an equation] : Calculating Taylor polynomial of an implicit function given by an equation, it does not work.

Basically the same approach as in this recent MSE thread.

The derivatives can be found iteratively using the constraint.

The set-up:

g[y_, z_] := (8 + 8 z - 4 z^3 - z^4 - y (1 + z)^4)/(1 + z);
poly[x_, y_, z_] := -8 x +
3 y + (-24 x + 12 y) z + (-24 x + 18 y) z^2 + (4 - 8 x +
12 y) z^3 + (3 + 3 y) z^4;
zconstraint = poly[x, y, z[x, y]];

First solve for z at x=0,y=0. We use the constraint only (and similar when solving for derivatives).

zsolns =
Solve[(zconstraint /. {x -> 0, y -> 0}) == 0, z[0, 0]]

(* Out= {{z[0, 0] -> -(4/3)}, {z[0, 0] -> 0}, {z[0, 0] ->
0}, {z[0, 0] -> 0}} *)

We'll work with the branch hitting -4/3.

zderivsolns =
Solve[{D[zconstraint, x] == 0,
D[zconstraint, y] == 0}, {D[z[x, y], x], D[z[x, y], y]}] /. {x ->
0, y -> 0} /. zsolns[]

(* Out= {{Derivative[1, 0][z][0, 0] -> 1/24,
Derivative[0, 1][z][0, 0] -> 1/192}} *)

Iterate to get the higher order derivatives needed for the series.

zderiv2solns =
Solve[{D[zconstraint, {x, 2}] == 0,
D[zconstraint, {x, 1}, {y, 1}] == 0,
D[zconstraint, {x, 2}, {y, 1}] == 0}, {D[z[x, y], {x, 2}],
D[z[x, y], {x, 1}, {y, 1}],
D[z[x, y], {x, 2}, {y, 1}]}] /. {x -> 0, y -> 0} /.
zsolns[] /. zderivsolns

(* Out= {{{Derivative[2, 0][z][0, 0] -> -(3/128),
Derivative[1, 1][z][0, 0] ->
-(11/3072), Derivative[2, 1][z][0, 0] -> 135/32768}}} *)

Now take a series for g, convert to a Taylor polynomial, and plug in values for z and the needed partial derivatives.

ser = Normal[Series[g[y, z[x, y]], {x, 0, 2}, {y, 0, 1}]];
ser /. zsolns[] /. zderivsolns /. zderivsolns /. zderiv2solns

(* Out= {{{{-(296/27) + x^2 (1/9 - (13 y)/1152) +
x (-(8/9) + y/36) - (2 y)/27}}}} *)
• Thanks for your help. But why do you " work with the branch hitting -4/3". In fact, in my problem there is the requirement, z->0 as x,y ->0. I have added this requirement in my question. If so, how can I expand the function? – Mark_Phys Mar 29 at 5:11
• I chose that one because it was simple and because I failed to notice that requirement in the post. Try it using the one where z goes to zero. Might work fine, if the multiplicity of the root does not cause trouble. – Daniel Lichtblau Mar 29 at 15:13
• Disregard, I see your point. It is a Puiseux series at z=0. – Daniel Lichtblau Mar 29 at 15:25