# Multivariable Taylor expansion does not work as expected

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:

Normal[Series[f[x,y],{x,0,2},{y,0,2}]]


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$:

Normal[Series[f[x,y],{x,0,1},{y,0,1}]]


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?

-

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.

-
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

-
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


-

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


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.

-

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"]

-
Not sure what you're getting at; Sum[] is perfectly capable of localizing indices… – J. M. May 6 '15 at 6:54