Calculating Taylor polynomial of an implicit function given by an equation

Question

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 1

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.

Answer 2

Let's try again. Start by defining the implicit equation

eqn = Sin[x y z[x, y]] + x + y + z[x, y] == 0

Let us first take derivatives of the condition F[x,y,z[x,y]]==0

n = 4;
eqns = Union[Flatten[Table[D[eqn, {x, i}, {y, j}], {i, 0, n}, {j, 0, n}]]] /. 
Thread[{x, y} -> 0];

and define a vector corresponding to the unknown derivatives

var = Union[Flatten[Table[
  D[z[x, y], {x, i}, {y, j}], {i, 0, n}, {j, 0, n}]]] /. 
Thread[{x, y} -> 0];

Now we solve for them (since eqn is always satisfied, all its derivatives should be zero)

sol = Solve[eqns, var][[1]]

and write the series solution to the implicit equation with these solutions:

Normal@Series[z[x, y] /. Thread[{ x, y} -> ϵ {dx, dy}],
{ϵ, 0, 4}] /. sol/. ϵ-> 1

(* dx dy (dx+dy)- (dx+dy) *)

which corresponds to the Taylor expansion of z[x,y] near zero satisfying eqn.

Does this answer your question?

EDIT

if you push this to 12th order you get

(* dx^3 dy^3 (dx+dy)-dx^2 dy^2 (dx+dy)-1/6 dx^3 dy^3 (dx^3+9 dx^2 dy+9 dx dy^2+dy^3)+1/3 dx^4 dy^4 (2 dx^3+9 dx^2 dy+9 dx dy^2+2 dy^3)+dx dy (dx+dy)-dx-dy *)

if you look for the formal solution for an arbitrary F[x,y,z] expanded around {a,b} it reads to first order (replacing in the above eqn=F[x,y,z[x,y]]==0 and expanding around {x,y}=={a,b})

 z(a,b)-(dx F^(1,0,0)(a,b,z(a,b))+dy F^(0,1,0)(a,b,z(a,b)))/F^(0,0,1)(a,b,z(a,b))

and to second order

 Normal@Simplify@(Series[
   z[x, y] /. Thread[{ x, y} -> {a, b} + ϵ {dx, dy}],
   {ϵ, 0, 2}] /. sol) /. ϵ -> 1 /. 
 Derivative[i_, j_, k_][F][__] :> Subscript[F, i, j, k]

to third order we have

Answer 3

I think you can use AsymptoticSolve for this. As in Jens' answer, I will scale both x and y by a scale factor ϵ:

AsymptoticSolve[Sin[x y z] + x + y + z == 0 /. {x->ϵ x, y->ϵ y}, {z, 0}, {ϵ, 0, 3}]

{{z -> (-x - y) ϵ + x y (x + y) ϵ^3}}

Setting ϵ to 1 produces the same answer as Jens'.

Answer 4

Calculating Taylor polynomial of an implicit function given by an equation can be done in these three simple steps:

1. we use Series to prepare the series of the function
2. we specify the function value at one point
3. we use SolveAlways to find the coefficients in the series so that the equation holds for all independent variables.

Here is the code for the example from the question:

z = Series[f[x, y], {x, 0, 2}, {y, 0, 2}];
f[0, 0] = 0;
z /. SolveAlways[Sin[x*y*z] + x + y + z == 0, {x, y} ]

Pavel