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?

• You're not using Mathematica notation. Have you used Mathematica before? Commented Apr 14, 2013 at 15:19
• Yes, I did. I didn't think Mathematica notation was necessary here, given I don't have any code to show at the moment. Commented Apr 14, 2013 at 15:22
• Try the function Series[ ]. Commented Apr 14, 2013 at 15:55
• @chris , this code doesn't specify the constraint that yields z[x,y], and probably as a consequence, it doesn't output the numeric value of all the partials of z in the series.. Any suggestions? Commented Apr 14, 2013 at 16:30
• No. It can be shown using implicit function theorem that the equation Sin[x y z] + x + y + z==0 (along with some other conditions) implies that there exists a rectangle centered at (0,0,0) in which z is defined as a function of x,y. However this function remains implicit. Using the same theorem it is possible to calculate partials of z at (0,0), and this allows to construct a Taylor polynomial of z[x,y] in that rectangle. The Taylor of degree 2 came out -y-x by manual calculation. I want to check if that is true. Commented Apr 14, 2013 at 16:44

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.

• well the good news is our solutions agree :-) Commented Apr 14, 2013 at 19:02
• Many thanks to both of you! The good news is the output agrees with what I calculated. Now what's left is to learn and understand your code :) Commented Apr 15, 2013 at 10:26
• Do you have any recommendations on a good concentrated learning resource for Multivariate Calculus programming in Mathematica? Commented Apr 15, 2013 at 10:37

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}]]] /.


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}]]] /.


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.

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

• I don't know Mathematica well enough to understand all that you did. It doesn't seem like what I meant, though. Maybe you can explain and elaborate a bit more... Commented Apr 14, 2013 at 17:23
• After you watch the Gruffalo :P Commented Apr 14, 2013 at 17:24

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'.

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