Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

Thanks in advance, Dror.

share|improve this question
1  
You're not using Mathematica notation. Have you used Mathematica before? –  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
1  
@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
show 1 more comment

2 Answers 2

up vote 4 down vote accepted

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.

share|improve this answer
    
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
add comment

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]

image

to third order we have

image

share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.