# 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(xyz)+x+y+z=0.

Any ideas?

-
You're not using Mathematica notation. Have you used Mathematica before? – Dr. belisarius Apr 14 '13 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. – Dror Apr 14 '13 at 15:22
Try the function Series[ ]. – bill s Apr 14 '13 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? – Dror Apr 14 '13 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. – Dror Apr 14 '13 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 :-) – chris Apr 14 '13 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 :) – Dror Apr 15 '13 at 10:26
Do you have any recommendations on a good concentrated learning resource for Multivariate Calculus programming in Mathematica? – Dror Apr 15 '13 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}]]] /.
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.

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... – Dror Apr 14 '13 at 17:23
After you watch the Gruffalo :P – Dror Apr 14 '13 at 17:24