Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

The basic multivariable Taylor expansion formula around a point is as follows:

$$ f(\mathbf r + \mathbf a) = f(\mathbf r) + (\mathbf a \cdot \nabla )f(\mathbf r) + \frac{1}{2!}(\mathbf a \cdot \nabla)^2 f(\mathbf r) + \cdots \tag{1}$$

In Mathematica, as far as I know, there is only one function, Series that deals with Taylor expansion. And this function surprisingly doesn't expand functions in the way the above multivariable Taylor expansion formula does. What I mean is that the function Series doesn't produce a Taylor series truncated at the right order.

For example, if I want to expand $f(x,y)$ around $(0,0)$ to order $2$, I think I should evaluate the following Mathematica expression:


But the result also gives order $3$ and order $4$ terms. Of course, I can write the expression in the following way to get a series truncated at order $2$:


but in this way I lose terms like $x^2$ and $y^2$, so it is still not right.

The formula $(1)$ gives each order in each term, so if the function Series would expand a function in the way formula $(1)$ does, there will be no problem.

I am disappointed that the Mathematica developers designed Series as they did. Does anyone know how to work around this problem?

share|improve this question
up vote 34 down vote accepted

It's true that the multivariable version of Series can't be used for your purpose, but it's still pretty straightforward to get the desired order by introducing a dummy variable t as follows:

Normal[Series[f[(x - x0) t + x0, (y - y0) t + y0], {t, 0, 2}]] /. t -> 1

$(x-\text{x0}) (y-\text{y0}) f^{(1,1)}(\text{x0},\text{y0})+\frac{1}{2} (x-\text{x0})^2 f^{(2,0)}(\text{x0},\text{y0})+(x-\text{x0}) f^{(1,0)}(\text{x0},\text{y0})+(y-\text{y0}) f^{(0,1)}(\text{x0},\text{y0})+\frac{1}{2} (y-\text{y0})^2 f^{(0,2)}(\text{x0},\text{y0})+f(\text{x0},\text{y 0})$

The expansion is done only with respect to t which is then set to 1 at the end. This guarantees that you'll get exactly the terms up to the total order (2 in this example) that you specify.

share|improve this answer
Thank you very much!! Your solution is elegant. simple and general. – matheorem Nov 22 '12 at 8:03

Here's an attempt:

multiTaylor[f_, {vars_?VectorQ, pt_?VectorQ, n_Integer?NonNegative}] :=
                   Sum[Nest[(vars - pt).# &, (D[f, {vars, \[FormalK]}]
                       /. Thread[vars -> pt]), \[FormalK]]/\[FormalK]!,
                       {\[FormalK], 0, n}, Method -> "Procedural"]

multiTaylor[f[x, y], {{x, y}, {0, 0}, 2}]
   f[0, 0] + y*Derivative[0, 1][f][0, 0] + x*Derivative[1, 0][f][0, 0] +
   (y*(y*Derivative[0, 2][f][0, 0] + x*Derivative[1, 1][f][0, 0]) +
    x*(y*Derivative[1, 1][f][0, 0] + x*Derivative[2, 0][f][0, 0]))/2
share|improve this answer
I admire your mathematica skill. Your answer is right. But I think the following Jens's solution is better, it is simple and more enlightening. What do you think? Thank you all the same!! – matheorem Nov 22 '12 at 8:21
No problem; you should always pick the answer that is most helpful to you and give that answer the acceptance, after all. I wanted to write something a bit more general (i.e. something that can also handle three, four... variables), so mine is a little bit complicated. – J. M. Nov 22 '12 at 8:24
Thanks for the function – it's very useful to me. One slight problem I just noticed occurs when f contains the variable k, which will be replaced by some integer after the expansion. To improve you function, you might want to wrap it in Module[{k}, (* function code *)] to protect the local variable. – David Zwicker Jul 30 '15 at 20:07
@David, you're right; I'll incorporate the fix later. Thank you! – J. M. Jul 30 '15 at 22:02

Another possibility is to strip the "too high" terms with a rule:

ser = Series[f[x, y], {x, x0, 2}, {y, y0, 2}];
Normal[ser] /.
 Derivative[m__][f][args__] /; Plus[m] > 2 :> 0

Mathematica graphics

share|improve this answer

Jens answered this nicely. But my preferred way of thinking about this question is via re-scaling (essentially using perturbation theory) instead of dummy variables.

For example, to generate the series expansion, re-scale all variables by s and expand by series coercion:

(f[x,y] /. Thread[{x,y} -> s {x,y}]) + O[s]^3

TraditionalForm output

This form is generally useful—and although one can use Normal, often working with the series is what you really end up wanting to do.

share|improve this answer

If f is dependent on a variable k, the above definition leads to a conflict. It is better to rename k in a less common name such as var.

multiTaylor[f_, {vars_?VectorQ, pt_?VectorQ, n_Integer?NonNegative}] :=
                   Sum[Nest[(vars - pt).# &, (D[f, {vars, var}] /. Thread[vars -> pt]),var]/var!,
                       {var, 0, n}, Method -> "Procedural"]
share|improve this answer
Not sure what you're getting at; Sum[] is perfectly capable of localizing indices… – J. M. May 6 '15 at 6:54

Your Answer


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.