series expansion - Calculating Taylor polynomial of an implicit function given by an equation

I'd like to write a procedure that will take

1. an equation: F(x,y,z) = 0


2. chosen variable: x


3. a point: (a,b)


4. degree: n

And the output, when exists, should be the Taylor polynomial of degree n of x as an implicit function of y,z given by F(x,y,z) = 0, around (a,b).

For example, calculate Taylor polynomial of degree 2 around (0,0) of z(x,y) , given by Sin[x y z] + x + y + z == 0.

Any ideas?

Answer

To be definite about what the goal is, I'm assuming you want the following result to appear:

Series[f[ x, y, z ] /. 
   z -> solution /. {x -> ϵ x, y -> ϵ y}, {ϵ, 0, 3}]

$O[\epsilon^4]$

This means the constraint function is zero to the desired order as a function of the variables x and y.

Here is a way to get this result:

f[x_, y_, z_] := Sin[x y z] + x + y + z


n = 3;

solution = Normal[
   Simplify@Series[
     Simplify[Normal[
        InverseSeries[
         Series[
          Normal[
            Series[
             f[ϵ x, ϵ y, ϵ z], {ϵ,

               0, n}]] /. ϵ -> 1, {z, 0, n}]]] /. {z -> 0, 
        x -> ϵ x, y -> ϵ y}], {ϵ, 0, 
      n}]] /. ϵ -> 1

(* ==> -x - y + x y (x + y) *)

In all the expansions I keep track of powers using ϵ, which is set to 1 at the end (see related answer here). The important step is to single out z as an expansion variable in f for which I then construct the inverse series and set it to zero (that's the step with z -> 0 where z actually stands for f because the series has been inverted). The last step is to again construct a series so that I get the powers of x and y nicely arranged.

With the resulting solution, you can check that the first equation that defined the problem is indeed satisfied.